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Abstract We explore self-similar hydrodynamic evolu- 
tion of central voids embedded in an isothermal gas of 
spherical symmetry under the self-gravity. More specifi- 
cally, we study voids expanding at constant radial speeds 
in an isothermal gas and construct all types of possible 
void solutions without or with shocks in surrounding en- 
velopes. We examine properties of void boundaries and 
outer envelopes. Voids without shocks are all bounded 
by overdense shells and either inflows or outflows in the 
outer envelope may occur. These solutions, referred to as 
type X void solutions, are further divided into subtypes 
Xi and Xii according to their characteristic behaviours 
across the sonic critical line (SCL). Void solutions with 
shocks in envelopes are referred to as type Z voids and 
can have both dense and quasi-smooth edges. Asymp- 
totically, outflows, breezes, inflows, accretions and static 
outer envelopes may all surround such type Z voids. 
Both cases of constant and varying temperatures across 
isothermal shock fronts are analyzed; they are referred 
to as types Zi and Zu void shock solutions. We apply 
the 'phase net matching procedure' to construct various 
self-similar void solutions. We also present analysis on 
void generation mechanisms and describe several astro- 
physical applications. By including self-gravity, gas pres- 
sure and shocks, our isothermal self-similar void (ISSV) 
model is adaptable to various astrophysical systems such 
as planetary nebulae, hot bubbles and superbubbles in 
the interstellar medium as well as supernova remnants. 
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1 Introduction 

Voids may exist in various astrophysical cloud or nebula 
systems of completely different scales, such as planetary 
nebulae, supernova remnants, interstellar bubbles and 
superbubbles and so forth. In order to understand the 
large-scale dynamic evolution of such gaseous systems, 
we formulate and develop in this paper an isothermal 
self-similar dynamic model with spherical symmetry in- 
volving central expanding voids under the self-gravity. In 
our dynamic model framework, a void is approximately 
regarded as a massless region with a negligible gravity 
on the surrounding gas envelope. With such idealization 
and simplification in our hydrodynamic equations, a void 
is defined as a spherical space or volume containing noth- 
ing inside. In any realistic astrophysical systems, there 
are always materials inside voids such as stellar objects 
and stellar winds or outflows etc. In Section 4.2.3, we 
shall show that in astrophysical gas flow systems that 
our isothermal self-similar void (ISSV) solutions are gen- 
erally applicable, such as planetary nebulae, interstellar 
bubbles or superbubbles and so on. 

Observationally, early-type stars are reported to blow 
strong stellar winds towards the surrounding interstellar 
medium (ISM). Hydrodynamic studies on the interac- 
tion of stellar winds with surrounding gases have shown 
that a stellar wind will sweep up a dense circumsteller 
shell (e.g. Pikel'ner & Shchcglov 1968; Avedisova 1972; 
Dyson 1975; Falle 1975). Such swept-up density 'wall' 
surrounding a central star thus form interstellar bubbles 
of considerably low densities inside (e.g. Castor, McCray 
& Weaver 1975). For example, the Rosette Nebula (NGC 
2237, 2238, 2239, 2246) is a vast cloud of dusts and gases 
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spanning a size of ~ 100 light years. It has a spectac- 
ular appearance with a thick spherical shell of ionized 
gas and a cluster of luminous massive OB stars in the 
central region, whose strong stellar winds and intense ra- 
diation have cleared a 'hole' or 'cavity' around the cen- 
tre and given rise to a thick spherical shell of ionized 
gases (e.g. Mathews 1966; Borland, Montmerle & Doom 
1986). Weaver et al. (1977) outlined a dynamic theory 
to explain interstellar bubbles. They utilized equations 
of motion and continuity with spherical symmetry. They 
gave an adiabatic similarity solution, which is applica- 
ble at early times and also derived a similarity solution 
including the effect of thermal conduction between the 
hotter (e.g. T lO^K) interior and the colder shell of 
swept-up ISM. Their solution was also modified to in- 
clude effects of radiation losses. Weaver et al. (1977) 
did not consider the self-gravity of gas shell which can 
be dynamically important for such a large nebula, and 
therefore possible behaviours of their self-similar solu- 
tions were fairly limited. For example, the thickness of 
the gas shell was limited to ~ 0.14 times of the bubble 
radius. In our model formulation, we ignore the gravity 
of the central stellar wind region and of star(s) embed- 
ded therein. Thus, this central region is treated as a void 
and we explore the self-similar dynamic behaviours of 
surrounding gas shell and ISM involving both the self- 
gravity and thermal pressure. Our ISSV solutions reveal 
that the gas shell of a cloud can have many more types 
of dynamic behaviours. 

Planetary nebulae (PNe) represent an open realm 
of astrophysical applications of our ISSV model, espe- 
cially for those that appear grossly spherical (e.g. Abell 
39; see also Abell 1966)0 During the stellar evolution, 
planetary nebulae emerge during the transition from the 
asymptotic giant branch (AGB) phase where the star 
has a slow AGB dense wind to a central compact star 
(i.e. a hot white dwarf) where it blows a fast wind in 
the late stage of stellar evolution (e.g. Kwok, Purton & 
Fitzgerald 1978; Kwok & Volk 1985; Chevalier 1997a). 
The high temperature of the compact star can be a 
source of photoionizing radiation and may partially or 
completely photoionize the dense slower wind. Chevalier 
(1997a) presented an isothermal dynamical model for 
PNe and constructed spherically symmetric global hy- 
drodynamic solutions to describe the expansion of outer 
shocked shell with an inner contact discontinuity of wind 
moving at a constant speed. In Chevalier (1997a), grav- 
ity is ignored and the gas flow in the outer region can be 
either winds or breezes. In this paper, we regard the in- 
ner expansion region of fast wind as an effective void and 
use ISSV solutions with shocks to describe dense shocked 
wind and AGB wind expansion. One essential difference 
between our ISSV model and that of Chevaher (1997a) 
lays in the dynamic behaviour of the ISM. In Chevalier 
(1997a), shocked envelope keeps expanding with a van- 

^ Most planetary nebulae appear elliptical or bipolar in 
terms of the overall morphology. 



ishing terminal velocity or a finite terminal velocity at 
large radii. By including the self-gravity, our model can 
describe a planetary nebula expansion surrounded by an 
outgoing shock which further interacts with a static, out- 
going or even accreting ISM. In short, the gas self-gravity 
is important to regulate dynamic behaviours of a vast gas 
cloud. Quantitative calculations also show that the lack 
of gas self-gravity may lead to a considerable difference 
in the void behaviours (see Section 4.1). Likewise, our 
ISSV model provides more sensible results than those of 
Weaver et al. (1977). We also carefully examine the in- 
ner fast wind region and show that a inner reverse shock 
must exist and the shocked fast wind has a significant 
lower expansion velocity than the unshocked innermost 
fast wind. It is the shocked wind that sweeps up the AGB 
slow wind, not the innermost fast wind itself. This effect 
is not considered in Chevalier (1997a). We also compare 
ISSV model with Hubble observations on planetary neb- 
ula NGC 7662 and show that our ISSV solutions are 
capable of fitting gross features of PNe. 

Various aspects of self-similar gas dynamics have been 
investigated theoretically for a long time (e.g. Sedov 1959; 
Larson 1969a, 1969b; Penston 1969a, 1969b; Shu 1977; 
Hunter 1977, 1986; Landau & Lifshitz 1987; Tsai & Hsu 
1995; Chevaher 1997a; Shu et al 2002; Lou & Shen 2004; 
Bian & Lou 2005). Observations also show that gas mo- 
tions of this kind of patterns may be generic. Lou & 
Cao (2008) illustrated one general polytropic example of 
central void in a self-similar expansion as they explored 
self-similar dynamics of a relativistically hot gas with 
a polytropic index 4/3 (Goldreich & Weber 1980; Fill- 
more & Goldreich 1984). The conventional polytropic 
gas model of Hu & Lou (2008) considered expanding 
central voids embedded in "champagne flows" of star 
forming clouds and provided versatile solutions to de- 
scribe dynamic behaviours of "champagne flows" in H 
II regions (e.g. Alvarez et al. 2006). In this paper, we 
systematically explore isothermal central voids in self- 
similar expansion and present various forms of possible 
ISSV solutions. With gas self-gravity and pressure, our 
model represents a fairly general theory to describe the 
dynamic evolution of isothermal voids in astrophysical 
settings on various spatial and temporal scales. 

This paper is structured as follows. Section 1 is an in- 
troduction for background information, motivation and 
astrophysical voids on various scales. Section 2 presents 
the model formulation for isothermal self-similar hydro- 
dynamics, including the self-similar transformation, an- 
alytic asymptotic solutions and isothermal shock con- 
ditions. Section 3 explores all kinds of spherical ISSV 
solutions constructed by the phase diagram matching 
method with extensions of the so-called "phase net" . In 
Section 4, we demonstrate the importance of the gas self- 
gravity, propose the physics on void edge and then give 
several specific examples that the ISSV solutions are ap- 
plicable, especially in the contexts of PNe and interstel- 
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lar bubbles. Conclusions are summarized in Section 5. 
Technical details are contained in Appendices A and B. 



2 Hydrodynamic Model Formulation 

We recount basic nonlinear Eulcr hydrodynamic equa- 
tions in spherical polar coordinates (r, 6, ip) with self- 
gravity and isothermal pressure under the spherical sym- 
metry. 



2.1 Nonlinear Euler Hydrodynamic Equations 
The mass conservation equations simply read 



dM 
~dt 



dM 
or 



, 



dM 
dr 



(1) 



where u is the radial flow speed, M(r, t) is the enclosed 
mass within radius r at time t and p{r^ t) is the mass 
density. The differential form equivalent to continuity 
equation ([T]) is 



dp 

m 



J_d_ 

j,2 Qj, 



(r^pu) = 



For an isothermal gas, the radial momentum equation is 



du 
'dt 



du 



U-— = — 



dr 



dp 
p dr 



GM 



(3) 



where G = 6.67 x 10"^ dyn cm^ g~^ is the gravitational 
constant, a = {p/ pY^"^ — [k^T /m)'^^'^ is the isothermal 
sound speed and p is the gas pressure, fee = 1-38 x 
10~^^ erg K^^ is Boltzmann's constant, T is the con- 
stant gas temperature throughout and m is the mean 
particle mass. Meyer (1997) and Chevalier (1997a) ig- 
nored the gas self-gravity in the momentum equation. 

A simple dimensional analysis for equations (HI) — ([3]) 
gives an independent dimensionless similarity variable 



r/{at) 



(4) 



involving the isothermal sound speed a. The consistent 
similarity transformation is then 



p(r, t) = a{x)/{AT:Gt^) 
M{r, t) = a^tm{x)/G , 



u{r, t) — av{x) 



(5) 



where Oi{x), m(x), v{x) are the dimensionless reduced 
variables corresponding to mass density p(r,t), enclosed 
mass M{r,t) and radial flow speed u{r,t), respectively. 
These reduced variables depend only on x (Shu 1977; 
Hunter 1977; Whitworth & Summers 1985; Tsai & Hsu 
1995; Shu et al 2002; Shen & Lou 2004; Lou & Shen 
2004; Bian & Lou 2005). Meyer (1997) adopts a different 
self-similar transformation by writing p = p{x) /r^ . By 
equation ([4]) above, we know that in Meyer (1997), p{x) 



is exactly equal to x"^ a{x)a^ / {'^nG) here. So the simi- 
larity transformation of Meyer (1997) is equivalent to 
similarity transformation ([5]) here but without the self- 
_gravity. Further analysis will show that transformation 
([5|) can satisfy the void boundary expansion requirement 
automatically (see Section 4.2.1). 

With self-similar transformation (U and equa- 
tion ID) yields two ordinary differential equations (ODEs) 

, , dm dm , /r.\ 

m-\- [v — X)—— — , — — = X a . (6) 
dx dx 

The derivative term dm/dx can be eliminated from these 
two relations in equation ([6]) to give 



m(x) — x'^a{x — v) 



(7) 



A nonnegative mass m(x) corresponds to a; — t; > 0. In 
figure displays of —v{x) versus x profiles, the lower left 
portion to the line v — x = is thus unphysical. We refer 
to the line v — x = as the "Zero Mass Line" (ZML); 
inside of x — v, there should be no mass, corresponding 
to a void. 

Relation ([7|) plus transformation ([4]) and ([5]) lead to 
two coupled ODEs from equations ^ and ([3]), namely 



[ix-vf-1] 



V? 1] 



dv 
dx 



1 da 
a dx 



a(x — v) 



{x — v) 



(8) 



a (x — v) 

X 



(x^v) (9) 



(Shu 1977). ODEs ^ and ^ differ from eqs (3) and (4) 
of Chevalier (1997a) by including the self-gravity effect. 

For ODEs ^ and the singularity at {x -v)'^ ^1 
corresponds to two parallel straight lines in the diagram 
of —v{x) versus x, representing the isothermal sonic crit- 
ical hues (SCL) (e.g. Shu 1977; Whitworth & Summers 
1985; Tsai k Hsu 1995; Shu et al. 2002; Lou & Shen 2004; 
Bian & Lou 2005). Asx — = — l<Ois unphysical for 
a negative mass, we have the SCL characterized by 



V = X — I , 



= 2/a 



(10) 



Two important global analytic solutions of nonlinear ODEs 
([5]) and ^ are the static singular isothermal sphere (SIS; 
e.g. Shu 1977) 



v = , 



m — 2x 



(11) 



and non-relativistic Einstein-de Sitter expansion solution 
for an isothermal gas 



a 



2 

3 ' 



9" 



(12) 



(e.g. Whitworth & Summers 1985; Shu et al. 2002). 

Let denote the value of x at a sonic point on SCL. 
A Taylor series expansion of v(x) and a{x) in the vicinity 
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of shows that solutions crossing the SCL smoothly at 
X* have the form of either 



-V = {1 - a;,) + 
2 



1 

2/3 



Ij {x - x») + 
1 I (a; - x») + 



(13) 



or 



-V — {1 — a;,) {x — a;*) + 

-(x — x*) + 



2 2 



(e.g. Shu 1977; Whitworth & Summers 1985; Tsai & Hsu 
1995; Bian & Lou 2005; see Appendix A of Lou & Shen 
2004 for higher-order derivatives). Thus eigensolutions 
crossing the SCL smoothly are uniquely determined by 
the value of a;* . Eigensolutions of type 1 and 2 are defined 
by equations p3|) and (fT4)l . Physically, equations p3)) 
and HH) describe how the gas behaves as it flows from 
subsonic to supersonic regimes across the SCL in the 
local comoving framework. 

An important numerical solution to ODEs ([8]) and ([9]) 
is the Larson-Penston (LP) solution (Larson 1969a; Lar- 
son 1969b; Penston 1969). This solution have an asymp- 
totic behaviour ti — > 2a;/3 and a 1.67 as x 0~^. 
And the LP solution is also an eigensolution of type 2 
(equation [T^ and passes through the SCL smoothly at 
3/ — 2.33. 

When X — > -|-oo, either at large radii or at very early 
time, solutions to two nonlinear ODEs ([8]) and ([9]) have 
asymptotic behaviours 



v = V+ 



2- A V {A/6-l){A-2) + 2V^/3 



X X^ X 

A A{2 -A) (4 - A)VA 
_ _i I 



,,3 



X 



m = Ax ~ AV + 



2x4 3x5 

A{2-A) A{A~4.)V 



model, the gas has a slower outgoing radial velocity and 
the density decreases more rapidly than that of Cheva- 
lier (1997a), because when the self-gravity is included, 
the gas tends to accumulate around the centre. 

Counterpart solutions to equations (fTT |) — (fT5|) can also 
be generalized for conventional and general polytropic 
gases (Lou & Wang 2006; Wang & Lou 2007; Lou & Cao 
2008; Wang & Lou 2008; Hu & Lou 2008). 



(14) 2.2 Isothermal Shock Jump Conditions 



For an isothermal fluid, the heating and driving of a 
cloud or a progenitor star (such as a sudden explosion of 
a star) at t = will compress the surrounding gas and 
give rise to outgoing shocks (e.g. Tsai & Hsu 1995). In 
the isothermal approximation, the mass and momentum 
should be conserved across a shock front in the shock 
comoving framework of reference 



Pd{ud - Us) = Pu{Uu - Us) 



(16) 



alpd + PdUdiud - Us) = alpu + PuUuiuu - Us) , (17) 

where subscripts d and u denote the downstream and 
upstream sides of a shock, respectively (e.g. Courant & 
Friedricks 1976; Spitzer 1978; Dyson & WiUiams 1997; 
Shen & Lou 2004; Bian & Lou 2005). Physically, we have 
Us — UdXsd — 0-uXsu = Ts/t as the outgoing speed of a 
shock with being the shock radius. Conditions p8|) 
and (Uni) below in terms of self-similar variables are de- 
rived from conditions p6|) and ()17|) by using the reduced 
variables w(x), x and a(x) and the isothermal sound 
speed ratio r = ad/a„ = Xsu/xsd, 



ctd/au = [vu - Xsu)/[r{vd - Xsd)] , 



(18) 



2x 



6x2 



(15) 

where V and A are two integration constants (e.g. Whit- 
worth & Summers 1985; Lou & Shen 2004), referred to as 
the reduced velocity and mass parameters, respectively. 
The non-relativisitic Einstein-de Sitter expansion solu- 
tion does not follow this asymptotic solution (fT5|) at large 
X. Chevalier (1997a) also presented asymptotic solutions 
of v{x) and a(x) at large x. Case A = 1 va. our solution 
(|15p should correspond to asymptotic behaviours of v{x) 
and a(x) in Chevalier's model (1997a; see his equations 
5 and 6). However, the coefficient of x"^ term for a(x) 
in his model is 1 while in our model, it is 1/2; and co- 
efficient of x^^ term for v{x) in his model is 2 while in 
our model, it is 2 — A. These differences arise by drop- 
ping the gravity in Chevalier (1997a). Physically in our 



Vd-Xsd-T{Vu-Xsu) = {TVd-Vu){Vu-Xsu){vd-Xsd) ■ 

(19) 

Consequently, we have r = {Td/T^)'^/^ with T being the 
gas temperature. Physics requires Td > T„ leading to 
T > 1. For T — 1, conditions (flS)) and ([19]) reduce to 

ad/au = {vu - Xs)/ivd - Xs) , (20) 



{vu - Xs){vd - Xs) = I , (21) 

where Xs = x„ = Xd is the reduced shock location or 
speed (e.g. Tsai & Hsu 1995; Chevalier 1997a; Shu et al. 
2002; Shen & Lou 2004; Bian & Lou 2005). 



5 



3 Isothermal Self-Similar Voids 

Various similarity solutions can be constructed to de- 
scribe outflows (e.g. winds and breezes), inflows (e.g. 
contractions and accretions), static outer envelope and 
so forth. These solutions will be presented below in order, 
as they are useful to construct ISSV solutions. 



3.1 Several Relevant Self-Similar Solutions 

We first show some valuable similarity solutions in ref- 
erence to earlier results of Shu (1977) and Lou & Shen 
(2004). These solutions behave differently as a; ^ -l-oo 
for various combinations of parameters V and A in asymp- 
totic solution (fT5|) . 

3.1.1 CSWCP and EWCS Solutions of Shu (1977) 

Shu (1977) presents a class of solutions: collapse solutions 
without critical point (CSWCP) and expansion- wave col- 
lapse solution (EWCS) (see Fig.[T]). 

The CSWCP solutions (light dotted curves in Fig. 
[J) have asymptotic behaviours of V = and A > 2 
according to solution (fT^ at large x, and describe the 
central free-fall collapse of gas clouds with contracting 
outer envelopes of vanishing velocities at large radii. 

The EWCS solution (the heavy solid curve in Fig. [1]) 
is obtained with A ^ 2+ and V = in solution (fT5|) . 
This solution is tangent to the SCL at = 1 and has 
an outer static SIS envelope (solution [TT|l and a free-fall 
collapse towards the centre ( Shu 1977; Lou & Shen 2004; 
Bian & Lou 2005); the central collapsed region expands 
in a self-similar manner. 

3.1.2 Solutions smoothly crossing the SCL twice 

Lou & Shen (2004) studied isothermal similarity solu- 
tions and divided Class 1 similarity solutions, which fol- 
low free-fall behaviours as a; ^ O"*", into three subclasses 
according to their behaviours at large x. Class la similar- 
ity solutions have positive V at large x (see solution [T5|). 
which describe a cloud with an envelope expansion. Class 
la solutions are referred to as 'envelope expansion with 
core collapse' (EECC) solutions by Lou & Shen (2004). 
The case of V = corresponds to Class lb solutions and 
V < renders Class Ic solutions. The CSWCP solutions 
belong to Class Ic and the EWCS solution belongs to 
Class lb. 

Lou & Shen (2004) constructed four discrete solutions 
crossing the SCL twice smoothly (i.e. satisfying equa- 
tions [13] and [T4|) by applying the phase diagram match- 
ing scheme (Hunter 1977; Whitworth & Summers 1985; 
see subsection 3.4 below). These are the first four exam- 
ples among an infinite number of discrete solutions (see 
Lou & Wang 2007). 
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Fig. 1 Isothermal self-similar solutions of Shu (1977): —v{x) 
versus x profile of EWCS (heavy solid curve) and a sequence 
of five CSWCP solutions (dotted curves) with mass parame- 
ter A values marked along these five curves {A = 2.2 to 3.0). 
Four solutions smoothly crossing the SCL (dash-dotted line to 
the upper right) twice of Lou & Shen (2004): three light solid 
curves labelled by numerals '1', '2' and '3' are type 2-type 1 
solutions with their smaller crossing points x,(l) matching 
the type 2 derivative and the larger crossing points x,{2) fit- 
ting the type 1 derivative. The light dotted curve labelled by 
'type 2-type 2' is the unique type 2-type 2 solution with both 
crossing points following the type 2 derivative (see Table 1 
for relevant parameters). The dash-dotted line to the lower 
left is the ZML of a; - n = 0. 



Let x*(l) and x,(2) be the smaller (left) and larger 
(right) cross points for each of the four solutions. As they 
are all Class I solution, a;,(l) is less than 1 and the be- 
haviours near a;,(l) are determined by type 2 (as defined 
by equationll4[) to assure a negative derivative d{—v)/dx 
at a;*(l) along the SCL. As in Lou & Shen (2004), we 
name a solution that crosses the SCL twice smoothly 
as 'type2-typel' or 'type2-type2' solution, which corre- 
sponds the derivative types at a;*(l) and a;*(2). Three of 
the four solutions of Lou & Shen (2004) are 'type2-typel' 
solutions and we further name them as 'type2-typel-l', 
'type2-typel-2' and 'type2-typcl-3' (see Fig.[T]and Table 
1 for more details). 



3.2 Isothermal Self-Similar Void (ISSV) Solutions 

Relation ([7]) indicates that the ZML v — x = separates 
the solution space into the upper-right physical part and 
the lower-left unphysical part in a ~v{x) versus x pre- 
sentation. 

For a solution (with x — v > Q) touching the ZML 
at xq, then v{xq) — xq holds on and so does m(xo) = 
there. Given the definitions of m(a:) = GM{r, t)l{a?t) 
and X = r/{at), condition m(xo) = indicates a spheri- 
cal isothermal gas whose enclosed mass M{atxo,t) van- 
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Table 1 Relevant parameters are summarized here for the 
four solutions crossing the SCL twice shown in Figure [1] 
Three of them cross the SCL at Xt{l) with the second kind 
of derivative and at x,{2) with the first kind of derivative. In 
the 'Solution' column of this table, they are named as 'type2- 
typel-1', 'type2-typel-2' and 'type2-typel-3' corresponding 
to the curves labelled by '1', '2' and '3', respectively. The 
fourth solution crosses the SCL at a;,(l) and a;, (2) both with 
the second kind of derivative and is named as 'type2-type2' in 
the 'Solution' column, corresponding to the curve labelled by 
'type 2-type 2'. The first and third type2-typel solutions and 
type2-type2 solution all belong to EECC solution (or Class la 
solution), while the type2-typel-2 solution belongs to Class Ic 
solution. In Lou & Shen (2004), this type2-type2 EECC solu- 
tion passes the SCL at x,{l) ^ 0.632' and a;'*(2) « 1.349. We 
reproduce these results. However, we here adopt a;t(l) ~ 0.71 
and x,(2) ~ 1.287, which are calculated with a higher accu- 
racy. 



Solution 


Class 


x,(l) 


x,(2) 


V 


A 


type2-typel-l 


EECC(Ia) 


0.23 


1.65 


1.8 


5.1 


type2-typel-2 


Ic 


2.5858 X 10"" 


0.743 


-0.77 


1.21 


type2-typel-3 


EECC (la) 


6 X 10"'' 


1.1 


0.3 


2.4 


type2-type2 


EECC(Ia) 


0.71 


1.287 


1.5 


4.7 



ishes, and thus a central void expands at a constant ra- 
dial speed axo- 

The condition v{xq) = xq marks a void boundary in 
expansion. Physically, xq is the start point of a stream- 
line as v(xo) — Xq indicates matters flowing outwards at 
a velocity of the boundary expansion velocity. However, 
a problem would arise when the gas density just outside 
the void boundary is not zero. For an isothermal gas, the 
central vacuum cannot resist an inward pressure from 
the outer gas across the void boundary. We propose sev- 
eral mechanisms/scenarios to provide sufficient energy 
and pressure to generate voids and maintain them for a 
period of time without invalidating our 'vacuum approx- 
imation' for voids. In Section 4, we present explanations 
and quantitative calculations about these mechanisms. 

Mathematically, to construct spherical isothermal self- 
similar voids is to search for global solutions reaching the 
ZML at Xq > 0. If a solution touches the ZML at point 
Xq, then Xq should be the smallest point of the solution 
with gas. Solutions with a negative mass are unphysical. 
So if w = X holds on at xq, both v{x) and a{x) should 
be zero in the range < a; < xq . 

Before exploring isothermal self-similar void (ISSV) 
solutions of spherical symmetry, we note a few properties 
of such solutions. Nonlinear ODEs ^ and ^ give the 
following two first derivatives at xq 



dv 
dx 



= 



da 
dx 



= 



(22) 



^ The vacuum solution v — and a = satisfies dimen- 
sional equations (HJ — ([3]); it is not an apparent solution to 
reduced ODE ([8]) because a common factor a has been taken 
out before we arrive at ODE (IS}. With this physical under- 
standing, we still regard v — and a = as a solution to 
ODEs © and ©. 



So in the —v{x) versus x and a{x) versus x presenta- 
tions, the right tangents to these solution curves at xo 
are horizontal. 

For spherical isothermal self-similar dynamics, we can 
show that across a void boundary, a{x) must jump from 
to a nonzero value (see Appendix A). Voids must be 
bounded by relatively dense shells with a density jump. 
We propose that the height of such a density jump may 
indicate the energy to generate and maintain such a void. 
Energetic processes of short timescales include super- 
novae involving rebound shock waves or neutrino emis- 
sions and driving processes of long timescales include 
powerful stellar winds (see Section 4 for more details). 
In reality, we do not expect an absolute vacuum inside a 
void. Regions of significantly lower density in gas clouds 
are usually identified as voids. 

For X — > -|-oo, the physical requirement of finite mass 
density and flow velocity can be met for a{x) and v(x) 
by asymptotic solution (fT5|) . 

ISSV solutions need to cross the SCL in the —v{x) 
versus x profiles as they start at the ZML and tend to 
a horizontal line at large x with a constant V. Given 
conditions (fTS]) and (fT4|) . ISSV solutions can be divided 
into two subtypes: crossing the SCL smoothly without 
shocks, which will be referred to as type and crossing 
the SCL via shocks, which will be referred to as type Z 
and can be further subdivided into types Zi and Zjj as 
explained presently. 



3.3 Type X ISSV Solutions without Shocks 

As analyzed in Section 3.2, type X ISSV solutions cross 
the SCL smoothly. Let x^, denote the cross point on the 
SCL. We have v{x.^,) = — 1 and a(x*) = 2/x* by 
equation (fTO|l . Conditions (fTS]) and (fT4|) give the eigen 
solutions for the first derivatives v'{x) and a'{x) at x* 
as either type 1 (equation [T51) or type 2 (equation fH)). 
Given x* and the type of eigen-derivative at x*, all the 
necessary initial conditions [w(x*), a(x*), w'(x*), a'(x*)] 
are available for integrating nonlinear ODEs ([8]) and ([9|) 
in both directions. 

We use Xq, ao = a(xo), x* and the types of eigen- 
derivative at x■^, to construct type X ISSV solutions. 
While Xo and ao are the key parameters for the void 
expansion speed and the density of the shell around the 
void edge, we use x* to obtain type X ISSV solutions as 
x^ parameter can be readily varied to explore all type X 
ISSV solutions. 

3.3.1 Type X\ ISSV Solutions: Voids with Sharp Edge, 
Smooth Envelope and Type 1 Derivative on the SCL 

Type Xi solutions cross the SCL smoothly and follow the 
Type 1 derivative at x* on the SCL. By equation (|13p . 
the first derivative d{—v)/dx is positive for < x* < 1 
and negative for x* > 1. This allows the x^ of type Xi 
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Table 2 Parameters of several typical type Xi ISSV so- 
lutions. Eight solutions under four different conditions dis- 
cussed in Section 3.3.1 are tabulated here. 



Condition 




Xo 




ao 


V 


A 


I 


0.10 


0.0132 




1.11 X lO'' 


-3.51 


0.070 


I 


0.50 


0.0417 




1.92 X 10^ 


-1.57 


0.649 


I 


0.70 


0.0238 




1.35 X 10^ 


-0.91 


1.10 


II 


0.75 


7.5 X 10~ 


5 


1.5 X 10"' 


-0.75 


1.23 


II 


0.85 


7.5 X 10" 


4 


2.5 X 10^ 


-0.46 


1.50 


II 


0.95 


5.6 X 10" 


4 


6.9 X 10^ 


-0.15 


1.83 


III 


1.05 


9.2 X 10~ 


5 


1.5 X 10^ 


-f0.14 


2.17 


IV 


1.60 


0.0242 




1.5 X 10" 


+1.71 


4.77 




X 



solutions to run from to -|-oo when behaviours of these 
solutions in the inner regions are ignored temporarily. 
When the existence of xq is required for constructing 
isothermal central voids, the range of a;, along the SCL 
is then restricted. 

Given on the SCL and using equations and 
p^ . we obtain the initial condition v(x^). q;(x*), 

to integrate ODEs ([8|) and Q from in both 
directions. If an integration towards small x can touch 
the ZML at a > and an integration towards -l-oo 
exists, then a type Xi ISSV solution is constructed. 

We now list five important numerals: xi w 0.743, 
X2 = 1, cca w 1.1, X4 « 1.65 and xs = 3. They are the 
cross points at the SCL of Lou & Shen type2-typel-2, 
SIS (see equation 11), Lou & Shen type2-typel-3, Lou & 
Shen type2-typel-l, Einstein-de Sitter solution, respec- 
tively. All these solutions have type 1 derivatives near 
their cross points on the SCL. However, their behaviours 
differ substantially as x 0+. The three solutions of 
Lou & Shen approach central free-fall collapses with a 
constant reduced core mass mo as the central mass ac- 
cretion rate; while SIS and Einstein-de Sitter solution 
have vanishing velocities as x ^ 0"*". 

Numerical computations show that type Xi ISSV so- 
lutions exist when x, falls into four intervals out of six 
intervals along x > axis divided by the five numer- 
als above. The four intervals are < x, < xi « 0.743, 
xi ~ 0.743 < X* < X2 = 1, X2 = 1 < X, < X3 « 1.1 and 
X3 « 1.1 < X, < X4 « 1.65. No type Xi exist with x* in 
intervals X4 w 1.65 < x, < X5 = 3 and x* > X5 = 3, be- 
cause integrations from x, in these two intervals towards 
-|-0 or -|-oo must halt when they encounter the SCL again, 
respectively. The six regions mentioned above are named 
as conditions I, II, III, IV, V and VI, respectively. Fig- 
ure [2] illustrates several typical type Xi ISSV solutions 
with their x* in different regions. The relevant solution 
parameters are summarized in Table 2. 




10"' 10"' 10"^ 10"' 10° 



X 

Fig. 2 Typical type Xi ISSV solutions. Panel A in linear 
scales shows —v{x) versus x curves of several type Xi solu- 
tions. The sonic critical point Xt of each curve is noted along 
the a;— axis. The dash-dotted curve which crosses the SCL at 
Xt = 1.7 smoothly and encounters the line again at ~ 0.4 
shows a typical behaviour when X4 « 1.65 < a;, < X5 = 3. 
The inset shows the enlarged portions of these ISSV solution 
curves near a; ^ 0^. The dash-dotted lines in both panel 
A and inset are the ZML. Panel B shows the log a versus 
X curves of the same solutions and the curves in panel A, 
inset and panel B with the same line type (light and heavy 
solid, dotted and dash) correspond to the same type Xi so- 
lutions. They are distinguished by their a;* values. Each of 
these curves jumps from a zero value (left) to a nonzero value 
(right) at xq, indicating a void in the region of a; < a;o. 



3.3.2 Type Xn Void Solutions: Voids with Sharp Edges, 
Smooth Envelope and Type 2 Derivative at the SCL 

Type Xii solutions cross the SCL smoothly and have 
the type 2 derivative at x*. By condition HH), the first 
derivative d{—v)/dx is negative for x* > 0. Because 
Xq < X*, X* of type Xii solutions must be larger than 
1 to assure behaviours of solution (|15p at large x (note 
that type 2 derivative is used to obtain a free-fall be- 
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Table 3 Parameters of several typical type Xu ISSV solu- 
tions. Five solutions under two conditions in Section 3.3.1 are 
listed below. 



0.5 



Condition 




Xo 


ao 


V 


A 


ir 


1.288 


2.8 X 10"* 


4 X 10** 


1.5 


4.7 


ir 


1.49 


0.022 


9.1 X 10=^ 


4.5 


18 


IV' 


2.35 


0.304 


1.61 


3.30 


8.54 


IV' 


2.5 


0.852 


1.32 


3.46 


8.92 


IV' 


2.58 


1.007 


1.23 


9.48 


3.54 



haviour around the centre). Similar to the approach to 
investigate type Xi ISSV solutions, we now list four im- 



portant numerals: x'l — 1, 



1.287, x'. = 1.50 and 



« 2.33. Here, X2 ~ 1.287 is the right cross point on the 
SCL of the type2-type2 solution of Lou & Shen (2004) 
and X4 « 2.33 is the cross point of the LP solution on 
the SCL. These two solutions both follow behaviours of 
type 2 derivative near their cross points. However, their 
behaviours differ as x 0+. Lou & Shen type 2-type 
2 EECC solution has a central free-fall collapse with a 
constant reduced core mass toq for the central mass ac- 
cretion rate; while LP solution has vanishing velocity 
and mass as a; — > O"*". And Xg = 1.5 is the critical point 
where the second-order type 2 derivative diverges (see 
Appendix A of Lou & Shen 2004). 

The four points given above subdivide the x > 1 
portion of the x— axis into four intervals: [x[ = l,xi, ^ 



1.287], 



1.287,4 = 1-5], W3 = 1.5,4 ~ 2.34] and 



[4 ~ 2.34, -|-cxd) which are referred to as condition F, IF, 
III', and IV', respectively. Similar to subsection 3.3.1, we 
choose a a;* and integrate from a;* in both directions. Nu- 
merical calculations indicate that type Xu void solutions 
only exist when their falls under conditions IF or IV'. 

We show four typical type Xu ISSV solutions in Fig- 
ure [3] with relevant parameters summarized in Table 3. 

3.3.3 Interpretations for Type X ISSV Solutions 

We have explored all possible type X ISSV solutions. 
Now we offer interpretations for these solutions. In pre- 
ceding sections, we used parameter on the SCL as 
the key parameter to construct type X ISSV solutions. 
This method assures us not to miss any possible ISSV 
solutions. However, for ISSV solutions, xq and ao are the 
most direct parameters describing properties of central 
void expansions. 

The void edge a;o is the reduced radius in the self- 
similar description and axQ is the expansion speed of the 
void boundary. We adjust the a;* value at the SCL to 
search for void solutions and there exists a certain rela- 
tionship between a;* and xq as expected. Figure |4] shows 
these relationships for types Xi and Xu ISSV solutions. 

Under conditions I, II, HI and given a:*, the cor- 
responding a;o runs from to a maximum value and 
then back to 0. The maximum values of conditions I, 



Sonic Critical Line 



Panel A 




0.01 0.02 0.03 



'Zero Mass line 



1 .288 1 .49 
X 



2.352.5 



10 - 

log(Q^ _ 



-5L 

10-' 







Panel B 




N. Sonic Critical Line >- , 




< 1 .288 








< 1 .49 



X 



0.5 - 


-0.5 - 



2.35 



Sonic Critical Line- 



Panel C 



-2.5 



1.5 

X 



Fig. 3 Four typical type Xu ISSV solutions in —v{x) ver- 
sus X profile. Values of their x, are marked by the grid lines 
of a;— axis in panel A. The inner box presents details of two 
curves near a; — > 0+ in the same line type as corresponding 
curves in Panel A. Panel B and C give the four solutions in 
log a versus x profile. Panel B presents two solution curves of 
ir condition marked with x^ = 1.288 (heavy solid line) and 
X, — 1.49 (light solid line) respectively and panel C shows 
two solution curves of IV' condition marked with x, — 2.35 
(light dash line) and x, = 2.5 (heavy dash line) respectively. 
The density jumps at xo are obviously much sharper in con- 
dition IF than in condition IV' (compare the numbers along 
y— axes). 



II and HI are xq — 0.042, xq — 8 x 10^'^ and xq = 
0.9 X lO^"' which correspond to « 0.5, a;, « 0.89 
and a;» w 1.04, respectivly. For a;* under conditions IV 
or IF, the corresponding xq increases with a:* and the 
intervals of xq are [0+, 0.027] and [0+, 0.022], respec- 
tively. Under condition IV', a;o increases with a;» from 
[0+, -|-cxd) monotonically. However, when x* is under con- 
dition IV', the corresponding xq is usually very large (at 
least 100 times larger) compared to other type X ISSV 
solutions unless x, is near 2.34. Physically under condi- 
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• Type )^ void solutions 
■o Type )^| void solutions 



-10' • ^ — --^ ' ' ' 1 

0.5 0.743 1 1.288 1.65 2 2.34 3 

a; » 

Fig. 4 The relation between isothermal void boundary xo 
and the sonic point x, at SCL. The 'solid circle curve' and 
the 'open circle curve' are the log(xo) versus x* curves of type 
Xi and type Xu void solutions. The key numerals dividing the 
positive real a;— axis into several ranges are indicated by the 
dash vertical lines with numerals marked along or at a;— axis. 
In regions I, II, III, and IV, we have four branches for type Xi 
void solutions. In regions IF and IV', we have two branches 
for type Xu ISSV solutions. 



tion IV', isothermal voids usually expand at a relatively 
high speed. Based on the value and meaning of xo, type 
X ISSV solutions can be divided into two classes: rapid 
and slow ISSV expansion solutions. All type Xi ISSV so- 
lutions and type Xu ISSV solutions under condition IF 
belong to slow void expansion solutions. Type Xu ISSV 
solutions under condition IV' belong to rapid void ex- 
pansion solutions. 

Parameter ao represents the reduced density a at the 
void boundary. Figures[2]and[3]and Tables 2 and 3 clearly 
indicate that isothermal voids, described by type X ISSV 
solutions, are all surrounded by dense mass shells and 
the gas density around voids attenuates monotonically 
with increasing radius. So there are classes of voids that 
evolve with fairly sharp edges but without shocks. Never- 
theless, sharp edge around voids is not a general property 
of all void solutions. Expanding voids with shocks can be 
surrounded by quasi-smooth edges (never smooth edges 
as shown in Appendix A). In Section 4, we will show 
that ao is an important parameter which may reveal the 
mechanism that generates and sustains a void. Large ao 
requires very energetic mechanisms against the high in- 
ward pressure across the boundary. So type X voids may 
be difficult to form because all type X voids have dense 
boundaries. 

Figures [2] and [3] and Tables 2 and 3 clearly show that 
all type Xu void solutions and type Xj ISSV solutions 
under conditions III and IV describe isothermal voids 
surrounded by gas envelopes in expansion (i.e. veloc- 
ity parameter I^ at a; — > cx) are positive). Astrophysical 
void phenomena are usually coupled with outflows (i.e. 



winds). Our ISSV solutions indicate that rapidly expand- 
ing voids must be surrounded by outflows. 

Type Xi ISSV solutions under conditions I and II 
describe voids surrounded by contracting envelopes, al- 
though under these two conditions the voids expand very 
slowly (< 0.042a, see subsection 3.3.1) and are surrounded 
by very dense shells (see Table 2). 

Outflows and inflows are possible as indicated by type 
X ISSV solutions, but no static shell is found in type X 
ISSV solutions. In the following section, we show voids 
with shocks being surrounded by static envelopes. 

A clarification deems appropriate here that division 
of the a; > axis in subsection 3.3.1 is actually not pre- 
cise. In subsection 3.3.1, we divide (0, -l-oo) into six in- 
tervals by five points xi to x^ and xi, X3 and X4 are 
the cross points at the SCL of Lou & Shen type2-typel 
solution. We note that they are only the first three exam- 
ples of an infinite number of discrete solutions that cross 
the SCL smoothly twice via type 2 derivative first at a 
smaller x and then type 1 derivative at a larger x (Lou & 
Shen 2004). Our numerical computations show that the 
fourth type2-typel solution will pass the SCL smoothly 
at x*(l) « 4 X 10"® and x*(2) = 0.97, so there should 
be another regime of 0.97 < x < 1 = X2 inside condi- 
tion II. The first four right cross points of type2-typel 
solutions are 1.65, 0.743, 1.1 and 0.97. By inference, the 
cross points of the fifth and following type2-typel solu- 
tions will be narrowly located around x = 1. When the 
infinite number of type2-typel solutions are taken into 
account, there will be fine structures around x — 1 in 
subsection 3.3.1 and Figure ID However, the solution be- 
haviours of crossing the SCL smoothly under condition 
of fine structure are like those under conditions II, III 
and IV near x — 1. 



3.4 Type Z Voids: ISSV Solutions with Shocks 

Shock phenomena are common in various astrophysical 
flows, such as planetary nebulae, supernova remnants, 
and even galaxy clusters gas (e.g. Castor et al. 1975; 
McNamara et al. 2005). In this subsection, we present 
type Z ISSV solutions, namely, self-similar void solu- 
tions with shocks. Equations and p?)) are mass and 
momentum conservations. The isothcrmality is a strong 
energy requirement. In our isothermal model, an exam- 
ple of polytropic process, the energy process is simplified. 
This simplification gives qualitative or semi-quantitative 
description of the energy process for a shock wave. By 
introducing parameter r for the temperature difference 
after and before a shock, we can describe more classes of 
shocks (Bian & Lou 2005). 

The basic procedure to construct a spherical ISSV so- 
lution with shocks is as follows. Given {xq, — xq, ao) 
at the void boundary, we can integrate ODEs ^ and ^ 
outwards from xq ; in general, numerical solutions cannot 
pass through the SCL smoothly (if they do, they will be 
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referred to as type X ISSV solutions); however, with an 
outgoing shock, solutions can readily cross the SCL (e.g. 
Tsai & Hsu 1995; Shu et al. 2002; Shen & Lou 2004; Bian 
& Lou 2005); finally global {xq < x < +oo) solutions can 
be constructed by a combination of integration from xq 
to Xds, shock jump and integration from Xus to +oo, 
where Xds and Xus are defined in subsection 2.2 as the 
radial expanding velocity of a shock on the downstream 
and upstream sides, respectively. 

A typical ISSV solution with a shock has four degrees 
of freedom (DOF) within a sensible parameter regime. 
For example, we need independent input of xq, ag, Xds 
and r to determine an ISSV solution with a shock, while 
the degree of freedom for type X solutions is one (i.e. 
a;* plus the type of eigen-derivative crossing the SCL 
are enough to make a type X ISSV solution). When we 
consider the simplest condition that r = 1, the DOF of a 
void solution with shocks is three. So infinite void shock 
solutions exist. By fixing one or two parameters, we can 
enumerate all possible values of the other parameter to 
obtain all possible ISSV solutions. For example, by fixing 
velocity parameter V, we can adjust mass parameter A 
and Xus to explore all void solutions (see equation I15[) ; 
following this procedure, xq, ao and Xds are determined 
hy V, A and Xus- There is a considerable freedom to set 
up an ISSV solution with shocks. In astrophysical flows, 
we would like to learn the expansion speed of a void, the 
density surrounding a void and the radial speed of gas 
shell at large radii. We then choose one or two parameters 
in Xq, ao and V as given parameters to search for ISSV 
solutions by changing other parameters such as Xus ■ 

In Section 3.4.1, we will first consider the simple case: 
equi-temperature shock void (i.e. r — 1), and refer to as 
void solutions with equi-temperature shocks or type Zi 
ISSV solutions. Several type Zi voids with different be- 
haviours near void boundaries and outer envelopes (a 
static envelope, outflows and inflows) will be presented. 
Phase diagram matching method will be described and 
extended to the so-called 'phase net', with the visual con- 
venience to search for ISSV solutions with more DOF. 
Section 3.4.2 presents type Zu ISSV solutions: void so- 
lutions with two-soundspeed shocks (i.e. r > 1). 



3.4-1 Type Z\ Void Solutions: Voids 
with Same-Soundspeed Shocks 

For an equi-temperature shock, we have t — 1 and thus 
Xds = Xus — Xs (see Section 2.2). As already noted, the 
DOF of Zi void solutions is three. We can use {xo, ao, 
Xs} to construct type Zi ISSV solutions. 

We have freedom to set the condition {xo, ao) at 
the void edge to integrate outwards. Before reaching the 
SCL, we set an equi-temperature shock at a fairly arbi- 
trary Xs to cross the SCL. We then combine the inte- 
grations from Xo to xj and from x^ to -|~oo by a shock 
jump to form a type Zi ISSV solution with a shock. 



We emphasize that under type Zi condition, the in- 
sertion of an equi-temperature shock does assure that 
void solutions jump across the SCL. Physics requires that 
at every point (x, v{x), a{x)) of any void solution, there 
must be a; > v{x) for a positive mass. Equation (|2ip indi- 
cates that across the equi-temperature shock front, the 
product of two negative Vd — Xs and w„ — Xs makes 1. 
In our model, Xs — Vd < 1 so Xs — Vu must be larger 
than 1. So the downstream and upstream are separated 
by the SCL. However, in type Zu condition with r > 1, 
this special property may not always hold on. We shall 
require such a separation across the SCL as a necessary 
physical condition. 

Figure [5] shows several type Zi ISSV shock solutions 
with void expansion at half sound speed and 0.03 times 
sound speed. From this figure, we know that even the 
voids expand at the same speed, with different density 
ao near the void edge and the radial velocities of the 
shock wave (xs), they can have outer envelopes of various 
dynamic behaviours. From values of v{x) and a{x) at 
large x (say 10), we can estimate V and A. Different 
from type X ISSV solutions, some type Zj ISSV shock 
solutions have outer envelopes with a negative V (e.g. 
curves 1' and 2' in Fig. [5]), that is, contracting outer 
envelopes. Numerical calculations show that voids with 
contracting envelopes usually expand very slowly (voids 
1' and 2' in Fig. [5]expand at 0.03a). 

In addition to shocks and outer envelope dynamics, 
another important difference between type Zj and type 
X ISSV solutions lies in behaviours of shells surround- 
ing the void edge. From Section 3.3 (see Fig. [2]), we have 
known that the reduced mass density of type X void solu- 
tions must encounter a sharp jump and decrease mono- 
tonically with increasing x. However, Fig. [5] (curves 4 
and 2' with the y-axis in logarithmic scales) and Fig. [6] 
indicate that with shocks involved, the density of shells 
near void edges can increase with increasing x. Under 
these conditions, density jumps from a void to gas ma- 
terials around voids appear not to be very sharp. Voids 
described by solutions like curves in Figure [S] and curves 
4 and 2' in Figure [5] have such 'quasi-smooth' edges. 
These solutions can approximately describe a void with 
a smooth edge, whose outer shell gradually changes from 
vacuum in the void to gas materials, without sharp den- 
sity jump. Or, they describe a void with a quasi-smooth 
edge, whose outer shell gradually changes from vacuum 
in the void to gas materials, with a small density jump. 
Figure [H] clearly shows that the faster a shock moves rela- 
tive to void edge expansion, the higher density rises from 
void edge to shock. 

Type Z\ Voids with a SIS Envelope Void phenomena in 
astrophysics indicate an expanding void in the centre 
and static gas medium around it in the outer space. For 
example, a supernova explodes and ejects almost all its 
matter into space. If the shock explosion approximately 
starts from the central core of the progenitor star, the 
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Fig. 5 We show six type Zi void shoclc solutions. Four of 
them have same xo = 0.5 but different values of qq marked 
by numbers 1, 2, 3, 4. The other two have same xo = 0.03 but 
different values of qq marked by number l',2',3',4'. Panel 
A above presents —v{x) versus x profiles and panel B be- 
low shows log a versus x profiles. Values of the two free pa- 
rameters (qo, Xs) axe (50, 0.8), (5, 1.5), (1, 2.3), (0.1, 2.5) 
for curves 1,2,3,4, respectively; and (100, 0.7), (0.01, 2) for 
curves 1' and 2', respectively. The dotted curves in both pan- 
els are Sonic Critical Line. 
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Fig. 6 Several type Zj ISSV solutions with quasi-smooth 
edges in light solid curves with the same values of a;o = 0.5 
and ao = 10~^ but different shock speed Xs from 2.5 to 5.0 
every 0.5 in step. Panel A presents a{x) versus x profiles and 
panel B shows —v{x) versus x profiles. The dotted and dash 
lines in panel B are the SCL and the ZML, respectively. 



remnant of the supernova is then approximately spheri- 
cally symmetric and a void may be generated around the 
explosion centre (e.g. Lou & Cao 2008). If the gravity of 
the central compact object may be ignored, we then de- 
scribe this phenomenon as an expanding spherical void 
surrounded by a static outer envelope. The analysis of 
Section 3.3.3 indicates that all type X ISSV solutions 
cannot describe this kind of phenomena. However with 
rebound shocks, it is possible to construct a model for 
an expanding void surrounded by a static SIS envelope. 

Shu (1977) constructed the expansion- wave collapse 
solution (EWCS) to describe a static spherical gas with 
an expanding region collapsing towards the centre. In 
fact, EWCS outer envelope with a; > 1 is the outer part 
of a SIS solution (see equation [TT|) . We now construct 
several ISSV solutions with an outer SIS envelope. An 
outer SIS envelope has fixed two DOF of a type Zi ISSV 
solution, that is, = and A — 2, so there is only one 
DOF left. A simple method is to introduce a shock at a 
chosen point Xs of EWCS solution except {x = 1, v = 
0, a — 1) (we emphasize that only one point (x = 1, v = 
0, a = 1) of the EWCS solution is at the SCL and all 
the other points lay on the upper right to the SCL in 
the —v(x) versus x profile) and make the right part of 
EWCS solution the upstream of a shock, then we can 
obtain (vds, ctds) on the downstream side of a shock. If 
the integration from (xg, Vds, ocds) leftward touches the 
ZML at Xo, a type Zi ISSV solution with a static outer 
envelope is then constructed. 

We introduce the a — v phase diagram to deal with 
the relationship among the free parameters and search 
for eigensolutions of DDEs dH) and Hunter (1977) in- 
troduced this method to search for complete self-similar 
eigensolutions of two DDEs. Whitworth & Summer (1985) 
used this method to combine free-fall solutions and LP 
solutions in the centre with certain asymptotic solutions 
at large radii. Lou & Shen (2004) applies this method to 
search for eigensolutions of ODEs ([5]) and ^ which can 
cross the SCL twice smoothly. In the case of Lou & Shen 
(2004), the DOF is 0. So there is an infinite number of 
discrete eigensolutions. 

For type Zi ISSV shock solutions with a static SIS 
outer envelope, the DOF is one. We insert a shock at 
Xg in the SIS and then integrate inwards from Xs to a 
fixed meeting point xp- Adjusting the value of Xs will 
lead to a phase curve [v(xp)'^, a^xp)^] in a versus v 
phase diagram. Meanwhile, an outward integration from 
a chosen void boundary condition (xq, vq — xq, ao) 
reaches a phase point [v{xp)~, a{xp)~] at xp. Vary- 
ing values of xq or ao will lead to another phase curve 
[v{xp)~ , a{xp)~] in the a versus v phase diagram. We 
note that changing both values of xq and ao will result 
in a "phase net" (i.e. a two-dimensional mesh of phase 
curves) in a versus v diagram (see Fig. [7]). If such a 
"phase net" of [xq, ao) and the phase curve of Xs share 
common points (usually, there will be an infinite num- 
ber of common points continuously as such "phase net" 
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Table 4 Values of xq, ao and Xs for several type Zi void 
shock solutions with an outer static SIS envelope are sum- 
marized here 
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Fig. 7 The phase diagram of a versus n at a chosen meet- 
ing point xf = 0.3 for match of the SIS solution (e.g., Shu 
1977) and type Z void shock solutions. Each open circle sym- 
bol joint by a heavy solid curve denotes an integration from a 
chosen Xa (marked besides each open circle) towards xf < Xs. 
We impose the equi-temperature shock conditions at Xs from 
0.32 to 1.36 for the outer SIS solution and then integrate from 
Xs back to Xf to get [v{xf), a{xF)] as marked in the phase 
diagram. The change of naturally leads to a phase curve 
shown here by the heavy solid curve. We choose the range 
of Xa to be [0.32, 1.36] because Xs must be larger than xf 
and a larger Xs than 1.36 will give rise to solution curves en- 
countering the SCL at a; > xf ~ 0.3. Meanwhile, two "phase 
nets", made of the light solid curves and dotted curves and 
connected by a medium heavy solid curve, are actually gen- 
erated by integrations from chosen (xq, vq = xq, ao). Each 
solid curve in the two nets (including the connecting medium 
heavy solid curve) is an equal— a;o curve, that is, every point 
in the same solid curve corresponds to the same value of xq 
noted besides each equal— a:;o curve; and each dotted curve 
in the two nets is an equal— ao curve with the value of ao 
in boldface noted besides each curve. For the lower right 
net, the points in the net correspond to an initial condi- 
tion of {xo, ao) e {xo\xo e [0,0.29]} x {ao|ao £ [5,10]}; 
the value of xo is 0, 0.15, 0.2, 0.25, 0.29 respectively from 
the left equal- a;o solid curve to the right one and the value 
of ao is 5, 7, 9, 10 respectively from the bottom equal- 
ao dotted curve to the top one. For the upper left net, 
the points in the net correspond to an initial condition of 
(xo, ao) e {a;o|a;o G [0,0.1]} x {ao|aa G [500,800]}; the value 
of Xo is 0, 0.04, 0.06, 0.08, 0.10 respectively from the right 
equal-xo solid curve to the left one and the value of ao is 500, 
600, 700, 800 respectively from the top equal-ao dotted curve 
to the bottom one. The medium solid curve is an equal- a;o, 
too, with a;o = and ao = 10, 50, 100, 200, 350 and 500 from 
lower right to left. This medium solid curve shows the trend 
of phase curves as we increase ao value in large steps. 
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is two dimensional), type Zi ISSV shock solutions with 
a static SIS envelope can be constructed. 

Figure [7] presents the phase diagram at a meeting 
point Xp = 0.3 to search for type Zi ISSV shock solu- 
tions with an outer SIS envelope. Note that part of the 
phase curve falls into the phase net, revealing that an in- 
finite number of type Zi ISSV shock solutions with outer 
SIS envelope indeed exist continuously. Shown by Figure 
[ll numerical results suggest that when shock position 
a;s > is less than 0.62 or larger than 1.335, there is at 
least one type Zi ISSV that can exist in the downstream 
side of a shock. However, if a shock expands at a radial 
velocity between 0.62a and 1.335a, it is impossible for a 
type Zi ISSV to exist inside a shock with a SIS envelope. 
Table 4 contains values of xq, ckq and Xg of some typi- 
cal type Zi ISSV shock solutions with a SIS envelope. 
Figure [8] is a phase diagram showing how xq and ao are 
evaluated with Xs changing to construct a type Zi ISSV 
shock solutions with an outer SIS envelope. 

All of Fig. [71 Fig. [8] and Table 4 clearly indicate type 
Zi ISSV shock solutions with a SIS envelope can be gen- 
erally divided into two classes according to Xg. Class I 
type Zi void shock solutions with an outer static SIS en- 
velope have Xs < 0.62 usually with a smaller value of xq 
and a higher value of ap. Class II solutions have Xg > 
1.335 usually with a larger value of xq and a medium 
value of ao- By a numerical exploration, the maximum 
of Xq of Class I solutions is ^ 0.09 for Xg = 0.26 and 
ao = 1.1 X 10'^; these voids expand at a low speed of 
< 0.1a; and the reduced density ao at the void bound- 
ary is usually > 10^, indicating a sharp edge density 
peak. Finally, we note that Class I voids involve shocks 
expanding at subsonic speeds. In this situation, the outer 
region x > Xg is not completely static. The SIS envelope 
only exists at x > 1, and the region between Xg and 1 is 
a collapse region (see two class I ISSV shock solutions 1 
and 2 in Fig. [9]). While in Class II ISSV shock solutions, 
shock expands supersonically and with xq usually rela- 
tively large. So the upstream side of a shock in Class II 
ISSV solutions is static (see two class II ISSV solutions 3 
and 4 shown in Fig. [9]). In rare situations, xq can be small 
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Fig. 8 The phase diagram of log ao versus log xo shows the 
relationship among Xs, xo and ao of type Zi void shock so- 
lutions with a static SIS envelope. The DOF of these three 
parameters is only one; i.e. given arbitrary one of the three, 
the other two parameters are determined in constructing a 
type Zi void shock solutions with a static SIS envelope. For 
any point in the two curves, xo and qq indicated by its x and 
y coordinates with Xs marked, correspond to a type Zi void 
shock solution. The upper left solid curve with its data points 
in asterisk symbol corresponds to the condition Xs < 0.62 and 
the lower dotted curve with its data points in open circle cor- 
responds to the condition Xs > 1.335. The first condition re- 
ferred to as Class I gives its largest xq = 0.09 when Xs = 0.26. 
Although Xo for the second condition referred to as Class II 
ranges along the entire real axis, it usually takes a relative 
large value; moreover, in Class II solutions, the reduced den- 
sity ao at the void boundary is in the order of unity. 



(e.g. Xs < 1.4) and ao neither large nor small, indicating 
that Class II voids have moderately sharp edges. 



Type Zi voids with expanding envelopes: breezes, winds, 
and outflows In our ISSV model, we use parameters V 
and A to characterize dynamic behaviours of envelopes. 
Equation (fTS)) indicates that > describes an expand- 
ing envelope at a finite velocity of Va, and a larger A 
corresponds to a denser envelope. For V = 0, the expan- 
sion velocity vanishes at large radii, corresponding to a 
breeze; a smaller A than 2 is required to make sure an 
outer envelope in breeze expansion. For V > 0^ the outer 
envelope is a wind with finite velocity at large radii. 

We apply the similar method to construct type Zi 
voids with expanding envelopes as we deal with type Zi 
voids with outer SIS envelopes. The difference between 
the two cases is that V and A are allowed to be different 
from V = and A — 2. In this subsection, we usually 
choose a meeting point xp between Xq and Xs- Then 
by varying xq and ao, we obtain a phase net composed 
by [v{xf)~, a(a;p')~]. Given V, we adopt A and inte- 
grate ODEs ([S]) and ^ from large x towards Xs- After 
setting a shock at Xs, we integrate ODEs towards xp- 
By varying A and Xg, we obtain a phase net composed 
by [v(xf)~^ , ct{xF)^]- The overlapped area of two phase 
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Fig. 9 Four typical type Zi void shock solutions with static 
SIS envelope. The two heavy solid curves in both panels are 
the Class I type Zi void solutions with a SIS envelope and a 
subsonic shock (i.e., Xa < 1), and the two light solid curves 
in both panels are the Class II type Zi void solutions with 
an outer SIS envelope and a supersonic shock (i.e., Xs > 1). 
Panel A presents —v{x) versus x profiles. Panel B presents a 
versus x profiles using a logarithmic scale along the y— axis. 
The key data (xo, ao, Xs) of these solutions are (0.09, 1.1 x 
10^ 0.26), (0.01, 9.5 X 10^ 0.62), (0.12, 7.5, 1.338), and 
(0.65, 4.1, 1.50) for curves 1, 2, 3, and 4, respectively. The 
dash curves in both panels are part of the EWCS solution. 

nets reveals the existence of type Zj ISSV with dynamic 
envelopes characterized by V and A. 

Figure [10] gives two examples of type Zi voids with 
breeze and wind. We emphasize that A > 2 is required 
in such ISSV solutions with an outer envelope wind and 
subsonic shock. Actually, the larger the velocity param- 
eter V is, the smaller the mass parameter A is needed. 
Large A is required to guarantee that the upstream re- 
gion of a global solution is on the upper right part of the 
SCL when the shock moves subsonically (i.e. Xg < 1). 
Physically, if an ISSV is surrounded by a subsonic shock 
wave, the wind outside needs to be dense enough. 

Type Z\ voids with contracting outer envelopes: accre- 
tions and inflows We explore ISSV with contracting en- 
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Table 5 Parameters of s ever al type Zi and Zn ISSV shock 
solutions shown in Figure [TT1 
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Fig. 10 Three typical type Zi void shock solutions with dif- 
ferent outer envelopes. Panel A presents —v{x) versus x pro- 
files. Panel B presents a{x) versus x profiles in a logarithmic 
scale along the y— axis. The heavy solid curve 1 gives a type 
Zi void shock solution with a quite thin breeze outer enve- 
lope whose (a;o, ao, Xs, V, A) = (1.62, 0.046, 2.6, 0, 0.1). 
The heavy dash curve 2 gives a type Zi void shock solution 
with a wind outer envelope whose {xo, ao, Xa, V, A) = 
(0.25, 28.2, 0.90, 0.2, 2.5). The heavy dotted curve 3 gives a 
type Zi void shock solution with an accretion outer envelope 
whose {xo, ao, Xs, V, A) = (0.50, 18.8, 1.00, 0, 2.50). The 
monotonic dotted curves in both panels stand for the SCL. 



velopes, such as accretion envelopes. Type Xi voids un- 
der conditions I and II have contracting envelopes. These 
type Xi voids are all surrounded by very dense shells with 
density decreasing with increasing radius. With shocks 
involved, central voids can have envelopes of various prop- 
erties. Type Zi voids with contracting outer envelopes 
are also studied. To have a contracting envelope, the ve- 
locity parameter V should be negative or approach 0~. 
A negative V and a positive A describe an outer enve- 
lope inflowing at a velocity of aV from large radii. For 
V = and A > 2, an outer envelope has an inflow ve- 
locity vanishing at large radii (Lou & Shcn 2004; Bian 
& Lou 2005). Figure [To] gives examples of type Zi voids 
with accreting outer envelopes. 

3.4-2 Type Z\\ ISSV Solutions: Voids Surrounded by 
Two-Soundspeed Shocks in Envelopes 

In the previous section, we explored type Zi ISSV solu- 
tions featuring the equi-temperature shock. There, r = 1 
indicates the same sound speed a across a shock. The 
isothermal sound speed a can be expressed as 



where ks is Boltzmann's constant, fi is the mean atomic 
mass and Z is the ionization state. Z = Q corresponds to 
a neutral gas and Z = 1 corresponds to a fully ionized 
gas. In various astrophysical processes, shock waves in- 
crease the downstream temperature, or change the pro- 
portion of gas particles; moreover, the ionization state 
may change after a shock passage (e.g. champagne flows 
in Hii regions, Tsai & Hsu 1995; Shu et al. 2002; Bian 
& Lou 2005; Hu & Lou 2008). Such processes lead to 
two-soundspeed shock waves with t > 1 . In this section, 
we consider r > 1 for type Zu ISSV shock solutions 
with two-soundspeed shocks. Global Zu void solutions 
have temperature changes across shock waves while both 
downstream and upstream sides remain isothermal, sep- 
arately. With a range of r > 1, it is possible to fit our 
model to various astrophysical flows. 

For T > 1, the DOF of type Zu ISSV shock solutions 
is four, i.e. one more than that of type Zi ISSV shock 
solutions with t — 1. We will not present details to con- 
struct type Zu ISSV shock solutions as they differ from 
the corresponding type Zi ISSV shock solution only in 
the quantitative sense. General properties such as the be- 
haviours near the void boundary and the outer envelope 
of type Zu ISSV shock solutions remain similar to those 
of type Zi ISSV shock solutions. We present examples of 
typical type Zu ISSV shock solutions in Figure [TTl 



4 Astrophysical Applications 

4.1 The role of self-gravity in gas clouds 

Earlier papers attempted to build models for hot bubbles 
and planetary nebulae (e.g. Weaver et al. 1977; Chevalier 
1997a) without including the gas self-gravity. In reference 
to Chevalier (1997a) and without gravity, our nonlinear 
ODEs dSD and ^ would then become 



dv 
dx 



-{x — v) 



(24) 



1/2 



{Z+l)ki 



■T 



1/2 



(23) 



1] = 
a dx X 



(25) 
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Fig. 11 Six type Z ISSV solutions with shocks. Panel A 
presents —vix) versus x profiles and panel B presents aix) 
versus x profiles. The heavy solid curves labelled by their cor- 
responding r in both panels are type Z ISSV shock solutions 
with = 2 and r = 1, 1.1, 1.2 respectively. The r = 1 
case is a type Zi ISSV shock solution and the latter two with 
r > 1 are type Z\\ ISSV shock solutions. The heavy dash 
curves labelled by their corresponding r in both panels are 
solutions with reduced velocity = and A = \, 1.5, 2.5 
respectively. The one with A = 2.5 has an envelope accre- 
tion and the other two have envelope breezes. The relevant 
data of these six solutions are summarized in Table 5. Shock 
jumps of type Zu ISSV solutions do not appear vertical as 
those of type Zi ISSV shock solutions (e.g. Bian & Lou 2005) 
because of different sound speeds across a shock front and 
thus different scales of reduced radial x = r /{at). 

ODEs p4)l and ([25)1 allow both outgoing and inflowing 
outer envelopes around expanding voids □ Self-similar so- 
lutions of ODEs ([24]) and ([25|l cannot be matched via 
shocks with a static solution of uniform mass density p. 
For comparison, the inclusion of self-gravity can lead to 
a static SIS. In some circumstances, there may be no 
apparent problem when ODEs and ([^5]) are applied 
to describe planetary nebulae because AGB wind outer 
envelope may have finite velocity at large radii. How- 
ever, in the interstellar bubbles condition, a static ISM 
should exist outside the interaction region between stel- 
lar wind and ISM. In Weaver et al. (1977), a static uni- 
formly distributed ISM surrounding the central bubble 
was specifically considered. In our ISSV model, a single 



^ Actually, Chevalier (1997a) did not consider physically 
possible situations of contracting outer envelopes. 
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Fig. 12 A comparison of ISSV solutions from different 
ODEs. Panel A presents —11(2;) versus x profiles and panel 
B presents a(x) versus x profiles. The light solid curves is 
a solution of ODEs ((Sjl and ((9)1 with self-gravi ty. T he heavy 
dashed curves represent a solution of ODEs (|24|l and (|25|l 
without self-gravity. An equi-temperature shock with r = 1 
is introduced in both solutions at Xs = 2.7025. The dotted 
curves in both panels represent the SCL. Other relevant pa- 
rameters are contained in Table 6. 



solution with shock is able to give a global description for 
ISM shell and outer region around an interstellar bubble. 
Naturally, our dynamic model with self-gravity is more 
realistic and can indeed describe expanding voids around 
which static and flowing ISM solutions exist outside an 
expanding shock front (Figs. [51 and [TTl Tables 4 and 5). 

For gas dynamics, another problem for the absence 
of gravity is revealed by asymptotic solution p5)) . where 
the coefficient of term in a{x) differs by a factor of 
1/2 between models with and without self-gravity, and 
the expression for v{x) at large radii differs from the 
a;~^ term. These differences lead to different dynamic 
evolutions of voids (see Section 2.1 for details). 

Under certain circumstances, the subtle difference be- 
tween the shell behaviours with and without gas self- 
gravity may result in quite different shell profiles around 
void regions. We illustrate an example for such differ- 
ences in Fig. [T2] with relevant parameters for the two 
solutions therein being summarized in Table 6. Given 
the same asymptotic condition of y = 0.2 and A = 1 
at large radii (Chevalier 1997a shows only A = 1 case), 
the behaviours of such voids differ from each other sig- 
nificantly with and without self-gravity, although both 
of them fit asymptotic condition (fT5|) well. With grav- 
ity, the given boundary condition and shock wave lead 
to a thin, very dense shell with a sharp void edge while 
the same condition leads to a quasi-smooth edge void 
without gravity. In Section 4.2, we will show that these 
two types of voids reveal different processes to generate 
and maintain them. In certain situations, void models 
without gravity might be misleading. 



16 



Table 6 Parameters for two ISSV solutions shown in Fig. 1121 
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4.2 Formation of ISSV Edge 

In our model, the central void region is simply treated as 
a vacuum with no materials inside. We need to describe 
astrophysical mechanisms responsible for creating such 
voids and for their local evolution. On the right side of 
the void edge Xq , the gas density a{x) > for x > Xo, 
while a{x) = for x < xq. li no materials exist within 
the void edge, there would be no mechanism to confine 
the gas against the inward pressure force across the void 
edge. We offer two plausible astrophysical scenarios to 
generate and maintain such voids. 

4-. 2.1 Energetics and Pressure Balance 

If we allow a tenuous gas to exist within a 'void' region to 
counterbalance the pressure across the 'void' edge for a 
certain time t, transformation ([5]) gives the edge density 
Po as 



pait) = 7.16 X 10% 



10^ yr 
t 



(26) 



where ao = ck^xq) is the reduced mass density at ISSV 
edge and rup is the proton mass. For a gas temperature 
T, the isothermal sound speed a is 



2.87 X 10^ 



T 



W K 



1/2 



(27) 



where the mean particle mass is that of hydrogen atom. 
Then the gas pressure po just on the outer side of the 
ISSV edge is 



Poit) = poa^ = 0.99q;o 



10^ yr 



T 



W K 



dyne cm , 
(28) 



where po is the mass density at the ISSV edge. Here, we 
take the proton mass as the mean particle mass. Equa- 
tion gives a pressure scaling po (x t~^T governed by 
self-similar hydrodynamics. 

Within a 'void', we may consider that a stellar wind 
steadily blows gas outwards with a constant speed. Var- 
ious astrophysical systems can release energies at differ- 
ent epochs. For an early evolution, massive stars steadily 
blow strong winds into the surrounding ISM (e.g. Math- 
ews 1966; Dyson 1975; Falle 1975; Castor et al. 1975; 
Weaver et al. 1977). In the late stage of evolution, com- 
pact stars can also blow fast winds to drive the surround- 
ing gas (e.g. Chevalier 1997a; Chevalier 1997b). As a 



model simplification, we assume that a tenuous gas mov- 
ing outwards at constant speed Vw with temperature T^, 
(we refer to this as a central wind) may carve out a cen- 
tral 'void' and can provide a pressure against the pres- 
sure gradient across the 'void' boundary; suppose that 
this central wind begins to blow at time t = 0. Then 
after a time t, the radius r of the central wind front is at 
fw = Vwt. By the mass conservation in a spherically sym- 
metric flow, the mass density of the central wind front 
at time t is then 



M 



Pw, front 



(29) 



where M is the mass loss rate of the central wind. For 
a contact discontinuity the ISSV edge between the in- 
ner stellar wind front and outer slower wind, the plasma 
pressure Pw, front of this central fast wind is 



Pw , front 



M 



which can be estimated by 

4.16 X 10-^M /10km s"^ 

IO-6M0 yr-i ^ 



lO^yr 



lO^K 



(30) 



dyne cm 



(31) 



where Mq is the solar mass and the mean particle mass 
is that of hydrogen atom. We adopt the parameters based 
on estimates and numerical calculations (see Section 4.3). 
By expressions ([50)1 and ([?T|) . the plasma pressure at the 
central wind front also scales as Pw, front oc t~^Tw By a 
contact discontinuity between the central wind front and 
the ISSV edge with a pressure balance Pw, front = Po, our 
self-similar 'void' plus the steady central wind can sus- 
tain an ISSV evolution as long as the central stellar wind 
can be maintained. Across such a contact discontinuity, 
the densities and temperatures can be different on both 
sides. Equations (PT|) and also show that, while pres- 
sures balance across a contact discontinuity, the reduced 
density at ISSV edge ao is determined by the mass loss 
rate M and central wind radial velocity Vw Back into 
this steady central stellar wind of a tenuous plasma at a 
smaller radius, it is possible to develope a spherical re- 
verse shock. This would imply an even faster inner wind 
inside the reverse shock (i.e. closer to the central star); 
between the reverse shock and the contact discontinuity 
is a reverse shock heated downstream part of the cen- 
tral stellar wind. Physically, the downstream portion of 
the central stellar wind enclosed within the contact dis- 
continuity is expected to be denser, hotter and slower 
as compared with the upstream portion of the central 
stellar wind enclosed within the reverse shock. 

The above scenario may be also adapted to supernova 
explosions. At the onset of supernova explosions, the flux 
of energetic neutrinos generated by the core collapse of 
a massive progenitor star could be the main mechanism 
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to drive the explosion and outflows. This neutrino pres- 
sure, while different from the central wind plasma pres- 
sure discussed above, may be able to supply sufficient 
energy to drive rebound shocks in a dense medium that 
trigger supernova explosions (e.g. Chevalier 1997b; Janka 
et al. 2007; Lou & Wang 2006, 2007; Lou & Cao 2007; 
Arcones et al. 2008; Hu & Lou 2009). Under certain con- 
ditions, such neutrino pressure may even counterbalance 
the strong inward pressure force of an extremely dense 
gas and generate central 'voids'. We will apply this sce- 
nario and give examples in Section 4.3. 

4.2.2 Diffusion Effects 

For astrophysical void systems with timescales of suf- 
ficient energy supply or pressure support being shorter 
than their ages, there would be not enough outward pres- 
sure at void edges to balance the inward pressure of the 
gas shell surrounding voids after a certain time. In such 
situations, the gas shell surrounding a central void will 
inevitably diffuse into the void region across the void 
edge or boundary. This diffusion effect will affect be- 
haviours of void evolution especially in the environs of 
void edge and gradually smear the 'void boundary'. How- 
ever, because the gas shell has already gained a steady 
outward radial velocity before the central energy supply 
such as a stellar wind pressure support fails, the inertia 
of a dense gas shell will continue to maintain a shell ex- 
pansion for some time during which the gas that diffuses 
into the void region accounts for a small fraction of the 
entire shell and outer envelope. As a result, it is expected 
that our ISSV solutions remain almost valid globally to 
describe void shell behaviours even after a fairly long 
time of insufficient central support. 

We now estimate the gas diffusion effect quantita- 
tively. We assume that a void boundary expands at ax^ 
and the void is surrounded by a gas envelope whose den- 
sity profile follows a p{r) oc fall-off (asymptotic so- 
lution [T5)) . The gas shell expands at radial velocity ax 
for X > xq. The central energy supply mechanism has 
already maintained the void to a radius tq and then fails 
to resist the inward pressure across the void edge. We 
now estimate how many gas particles diffuse into void 
region in a time interval of At ~ ro/(axo), during which 
the void is supposed to expand to a radius 2ro. We erect 
a local Cartesian coordinate system in an arbitrary vol- 
ume element in gas shell with the cc— axis pointing ra- 
dially outwards. The Maxwellian velocity distribution of 
thermal particles gives the probability density of velocity 
^ = {vx,Vy,Vz) as 

p„(f;)ocexp — . (32) 

Define I as the mean free path of particles in the gas shell 
near void edge. If a gas particle at radius r can diffuse 
into radius f without collisions, its velocity is limited by 



(r -f VxAtY + {vyAt)"^ + {vzAtY < r'^ and its position is 
limited by r — f < I. We simplify the velocity limitation 
by slightly increasing the interval as (r -I- VxAt)"^ < 
and Vy + v^< (r/At)'^. We first set r = tq and integrate 
the ratio of particles that diffuse into radius ro during 
At to total gas shell particles within radius vq + I a,s 

Iro 47rr^p(r)(ir Jllr+v^At\\<roizv2+v^<(ro/Aty Pvd^v 

_ 1 - exp [-rg/(2a2z\t2)] 
" (27r)i/2^ 

X dr exp {-v'^/2)dv 

"''■0 "'-^-^0 
_ 1 - exp (-x§/2) rp 
(27r)i/2 I 

pl+l/ra I'-x 

X dx exp {-i^/2)di , (33) 

J 1 J —x — 2xq 

where both x and v are integral elements. We simply 
set Xo = 1 and present computational results in Table 
7. It is clear that even the inner energy supply fails to 
sustain inward pressure for a fairly long time, there are 
only very few particles that diffuse into the original void 
region, namely, a void remains quite empty. 

In the context of PNe, the particle mean free path 
I may be estimated for different species under various 
situations. For an example of PN to be discussed in Sec- 
tion 4.3.1, / = l/(n(T) = 3 X lO^o cm, where n « 5000 
cm^'^ is the proton (electron) number density in the H II 
region and a = 6.65 x 10~^^ cm^ is the electron cross sec- 
tion in Thomson scattering. One can also estimate cross 
section of coulomb interaction between two protons as 
~ 10~^^ cm^ and thus the mean free path for proton 
collisions is ^ 10^3cm. A PN void radius is ~ 5 x 10^^ 
cm. If no inner pressure is acting further as a void ex- 
pands to radius ^ 10^^ cm, gas particles that diffuse into 
r < 5 X 10^^ cm only take up ~ 6 x 10~^ of those in the 
gas shell. However, at the onset of a supernova explosion, 
the particle mean free path in the stellar interior is very 
small. If the density at a void edge is 1.2 x 10^ g cm^"^ 
(see Section 4.3.2), and the scattering cross section is es- 
timated as the iron atom cross section 4.7 x 10~^* cm^, 
the particle mean free path is only / ^ 10^^^ cm. Under 
this condition, particle that diffuse into a void region can 
account for 6.2% of those in the total gas shell. 

In short, while diffusion effect inevitably occurs when 
the inner pressure can no longer resist the inward gas 
pressure across a void edge, it usually only affects the 
gas behaviour near the void edge but does not alter 
the global dynamic evolution of gas shells and outer en- 
velopes over a long time. However, we note that for a 
very long-term evolution, there will be more and more 
particles re-entering the void region when no sufficient 
pressure is supplied, and eventually the diffusion effect 
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Table 7 Ratio of molecules that diffuse into radius ro during 
At to total gas shell molecules within radius ro + I under 
different ratio l/ro. We take xq = 1. 
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will result in significant changes of global dynamical be- 
haviours of voids and shells. These processes may happen 
after supernova explosions. When neutrinos are no longer 
produced and rebound shocks are not strong enough to 
drive outflows, a central void generated in the explosion 
will gradually shrink in the long-term evolution of super- 
nova remnants (SNRs). 

4-2.3 Applications of IS SV Model Solutions 

In formulating the basic model, we ignore the gravity of 
the central void region. By exploring the physics around 
the void boundary, a tenuous gas is unavoidable inside 
a 'void' region for either an energy supply mechanism 
leading to an effective pressure or that diffused across 
the void boundary. As long as the gravity associated with 
such a tenuous gas inside the 'void' is sufficiently weak, 
our ISSV model should remain valid for describing the 
large-scale dynamic evolution of void shells, shocks, and 
outer envelopes. 

In a planetary nebula or a supernova remnant, there 
is usually a solar mass compact star at the centre (e.g. 
Chevalier 1997a). For outgoing shells at a slow velocity 
of sound speed ~ 10 km s~^, the Parker-Bondi radius 
of a central star of IOMq is ^ 10^^ cm (i.e. ^ 10"^ 
light year or ~ 3 x 10^'* pc). The typical radius of a 
planetary nebula is about 10^* cm (see Section 4.3.1). 
Even the youngest known supernova remnant G1.9-f0.3 
in our Galaxy, estimated to be born 140 yr ago, has a 
radius of - 2 pc (e.g. Reynolds 2008; Green et al. 2008). 
Thus the central star only affects its nearby materials 
and has little impact on gaseous remnant shells. 

For a stellar wind bubble (e.g. Rosette Nebula), there 
usually are several, dozens or even thousands early-type 
stars blowing strong stellar winds in all directions. For 
example, the central 'void' region of Rosette Nebula con- 
tains the stellar cluster NGC 2244 of ~ 2000 stars (e.g. 
Wang et al. 2008). Conventional estimates show that 
the thick nebular shell has a much large mass of around 
10,000-16,000 solar masses (e.g. Menon 1962; Krymkin 
1978). For a sound speed of ~ 10 km s~^, the Parker- 
Bondi radius of a central object of 2OOOM0 is then ~ 0.08 
pc, which is again very small compared to the ~ 6 pc 
central void in Rosette Nebula (e.g. Tsivilev et al. 2002). 
For a typical interstellar bubble consdiered in Weaver et 
al. (1977), the total mass inside the bubble, say, inside 
the dense shell, is no more than 50 solar masses, which 
is significantly lower than that of the 2000 solar mass 
dense shell. 



We thus see that the dynamical evolution of flow sys- 
tems on scales of planetary nebulae, supernova remnants 
and interstellar bubbles are only affected very slightly by 
central stellar mass objects (e.g. early-type stars, white 
dwarfs, neutron stars etc). Based on this consideration, 
we regard the grossly spherical region inside the outer 
dense shells in those astrophysical systems as a void in 
our model formulation, ignore the void gravity, empha- 
size the shell self-gravity and invoke the ISSV solutions 
to describe their dynamic evolution. 

4.3 Astrophysical Applications 

Our ISSV model is adaptable to astrophysical flow sys- 
tems such as planetary nebulae, supernova explosions, 
supernova remnants, bubbles and hot bubbles on differ- 
ent scales. 

4-3.1 Planetary Nebulae 

In the late phase of stellar evolution, a star with a main- 
sequence mass < 8Mq makes a transition from an ex- 
tended, cool state where it blows a slow dense wind to a 
compact, hot stat^ where it blows a fast wind. The inter- 
action between the central fast wind with the outer slow 
wind results in a dense shell crowding the central region 
which appears as a planetary nebula (e.g. Kwok, Purton 
& Fitzgerald 1978). The hot compact white dwarf star at 
the centre is a source of photoionizing radiation to ion- 
ize the dense shell (e.g. Chevalier 1997a). When the fast 
wind catches up the slow dense wind, a forward shock 
and a reverse shock will emerge on outer and inner sides 
of a contact discontinuity, respectively. As shown in Sec- 
tion 4.2.3, the gravity of central white dwarf and its fast 
steady wind is negligible for the outer dense wind. Prac- 
tically, the region inside the contact discontinuity may be 
regarded approximately as a void. Meanwhile, the pho- 
toionizing flux is assumed to be capable of ionizing and 
heating the slow wind shell to a constant temperature 
(Chevalier 1997a), and the outer envelope, the cool AGB 
slow wind that is little affected by the central wind and 
radiation, can be also regarded approximately as isother- 
mal. The constant temperatures of dense photoionized 
shell and outer envelope are usually different from each 
other which can be well characterized by the isothermal 
sound speed ratio t (see Section 2.2). Thus the dynamic 
evolution of dense shell and outer AGB wind envelope 
separated by a forward shock is described by a type Z 
ISSV solution. 

Within the contact discontinuity spherical surface, 
there is a steady downstream wind blowing outside the 
reverse shock front (e.g. Chevalier & Imamura 1983). 
Consistent with Section 4.2.1, we define as the ra- 
dius where the downstream wind front reaches and as 

■* A hottest white dwarf (KPD 0005 5106) detected recently 
has a temperature of ~ 2 x 10^ K (Werner et al. 2008) . 
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the radius of reverse shock. Within the radial range of 
< r < r^, i.e. in the downstream region of the reverse 
shock, we have a wind mass density 



Pwir, t) 



M 



(34) 



We define a^, and T^^d{u) a-s the sound speed and gas 
temperature on the downstream (upstream) side of the 
reverse shock and ratio = awA/aw,u = {Twa/Tw^uY/'^ 
to characterize the reverse shock. For a reverse shock in 
the laboratory framework of reference given by shock 
conditions (1161) and (1171). we have in dimensional forms 
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where the outgoing speed of the reverse shock, 

Uu,rs is the upstream wind velocity, Pu{d).rs is the up- 
stream (downstream) mass density, respectively. The first 
and second expressions in equation (|35p are equivalent: 
the first expresses the upstream flow velocity in terms of 
the downstream parameters, while the second expresses 
the downstream flow velocity in terms of the upstream 
parameters. In solving the quadratic equation, we have 
chosen the physical solution, while the unphysical one 
is abandoned. Normally, in the downstream region < 
r < rw, the plasma is shock heated by the central faster 
wind with r^, > 1. In the regime of an isothermal shock 
for effective plasma heating, we take — 1 (i.e. aw,d = 
o.w,u — a-iii) for a stationary reverse shock in the labora- 
tory framework of reference, shock conditions ([35]) and 
(l34l) reduce to 
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In this situation, the reverse shock remains stationary in 
space and this may shed light on the situation that an 
inner fast wind encounters an outer dense shell of slow 
speed. While a reverse shock always moves inwards rel- 
ative to both upstream and downstream winds, either 



outgoing or incoming reverse shocks are physically al- 
lowed in the inner wind zone in the laboratory frame- 
work of reference. In the former situation, the reverse 
shock surface, contact discontinuity surface and forward 
shock surface all expand outwards steadily, with increas- 
ing travelling speeds, respectively. In the latter situa- 
tion, the downstream wind zone, namely, the shocked 
hot fast wind zone expands both outwards and inwards, 
and eventually fills almost the entire spherical volume 
within the dense shell. The reverse shock here plays the 
key role to heat the gas confined within a planetary neb- 
ula to a high temperature, which is thus referred to as a 
hot bubble. In all situations, between the reverse shock 
and the contact discontinuity, the downstream wind has 
a constant speed, supplying a wind plasma pressure to 
counterbalance the inward pressure force across the con- 
tact discontinuity. 

Here, type Z ISSV solutions arc utilized to describe 
the self-similar dynamic evolution of gas shell outside 
the outgoing contact discontinuity. An outgoing forward 
shock propagates in the gas shell after the central fast 
wind hits the outer dense shell. According to properties 
of type Z ISSV solutions, there may be outflows (i.e. 
winds and breezes), static ISM or inflows (i.e. accretion 
flows and contractions) in the region outside the forward 
shock. The spatial region between the contact disconti- 
nuity and the forward shock is the downstream side of 
the forward shock. 

We now provide our quantitative model PN estimates 
for a comparison. Guerrero et al. (2004) probed the struc- 
ture and kinematics of the triple-shell planetary nebula 
NGC 7662 based on long-slit cchelle spectroscopic ob- 
servations and Hubble Space Telescope archival narrow- 
band images. They inferred that the nebula with a spa- 
tial size of ~ 4 X 10^* cm consists of a central cavity 
surrounded by two concentric shells, i.e. the inner and 
outer shells, and gave a number density Hp distribution 
(the dotted curve in panel A of Figure fT3)) . The temper- 
atures of the inner and outer shells were estimated as 

1.4 X 10"* K and ~ 1.1 x lO"' K, respectively. No in- 
formation about the inner fast wind is given in Guerrero 
et al. (2004). In our model consideration, the planetary 
nebula NGC 7662 may be described by our type Z ISSV 
model with a shocked inner fast wind. In our scenario 
for a PN, the central cavity in the model of Guerrero 
et al. (2004) should actually involve an inner fast wind 
region with a reverse shock. The inner and outer shells 
correspond to the downstream and upstream dense wind 
regions across a forward shock, respectively. Thus the 
inner boundary of the inner shell is the contact disconti- 
nuity in our model scenario. Physically, we suppose that 
the central star stops to blow a dense slow wind of ~ 10 
km s~^ about ~ 1000 years ago and the inner fast wind 
of 10^ K began to blow outwards from a white dwarf at 
'Uu,rs — 1500 km s~^ about ~ 600 years ago. When the 
inner fast wind hits the dense slow wind ~ 4 years after 
its initiation, a reverse shock and a forward shock are 
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Fig. 13 A type Z ISSV shock model with an inner fast wind 
to fit the planetary nebula NGC 7662. Panels A, B and C 
show our model results for proton number density np, radial 
flow velocity and temperature of NGC 7662 in solid curve, 
respectively. Numerals 1, 2 and 3 in Panel A mark reverse 
shock, contact discontinuity and forward shock surface, re- 
spectively. The dashed curve in Panel A is the estimate of 
proton number density by Guerrero et al. (2004). The in- 
ner fast wind blows at 1500 km s~^ inside the reverse shock, 
which is not shown in Panel B. 



generated on the two sides of the contact discontinuity. 
The reverse shock moves inwards at a speed ~ 10 km 
s""'^. The best density fit to the estimate of Guerrero et 
al. (2004) is shown in Figure [T51 In our model, the inner 
fast wind has a mass loss rate from the compact star as 
~ 2 X 1O~®M0 yr~^ (consistent with earlier estimates of 
Mellema 1994 and Perinotto et al. 2004) and the reverse 
shock is able to heat the downstream wind to a temper- 
ature of ~ 6.4 X 10^ K. The downstream wind of the 
reverse shock has an outward speed of ^ 25.9 km s~^, 
corresponding to a kinematic age (i.e. the time that a 
shocked fast wind at that velocity blows from the central 
point to its current position) of ^ 630 years. In Guerrero 
et al. (2004) , the kinematic age is estimated to be ~ 700 
years. In Guerrero et al. (2004), the inner shell density is 
~ 5 X 10'^ mp cm~^ and our model shows a density varia- 



tion from ~ 4 X lO"^ 



cm"'^ to ~ 7x lO"* 



with 



a comparable mean. The forward shock travels outwards 
at a speed ~ 43.0 km s~^, consistent with an average 
outward velocity of the inner shell at ~ 44 km s~^. The 
total mass of the inner shell is ~ 8.5 x 1O~'^M0, and the 
mass of the outer shell within a radius of 2.5 x 10^^ cm 
is ^ 0.036Mq, which are all consistent with estimates of 
Guerrero et al. (2004). However, Guerrero et al. inferred 
that the outer shell has an outward velocity of around 50 
km s~^ and a proton number density of ^ 3000 cm^'^. 



Our model estimates indicate that the outer shell has a 
proton number density from 3200 cm~^ at the immedi- 
ate upstream side of the forward shock to ~ 400 cm~'^ at 
2.5 X 10"'^* cm, and the outward velocity varies from ~ 26 
km at the upstream point of the forward shock to 
~ 20 km at 2.5 x 10^^ cm. And thus the dense slow 
wind mass loss rate is 0.68 x 10~^MQyr^^, which is con- 
sistent with earlier numerical simulations (e.g. Mellema 
1994; Perinotto et al. 2004), but is lower by one order 
of magnitude than ^ lO~^M0yr~^ estimated by Guer- 
rero et al. (2004). In summary, our ISSV model appears 
consistent with observations of the NGC 7662, and a 
combination of hydrodynamic model with optical and 
X-ray observations would be valuable to understand the 
structure and dynamic evolution of planetary nebulae. 

J1-.3.2 Supernova Explosions and Supernova Remnants 

At the onset of a type 11 supernova (or core-collapse su- 
pernova) for a massive progenitor, extremely energetic 
neutrinos are emitted by the neutronization process to 
form a 'neutrino sphere' that is deeply trapped by the 
nuclear-density core and may trigger a powerful rebound 
shock breaking through the heavy stellar envelope. At 
that moment, the central iron core density of a ~ 15Mq 
progenitor star can reach as high as ~ 7.3 x 10® g cm""^ 
and the core temperature could be higher than ~ 7.1 x 
10® K. The density of the silicon layer is ~ 4.8 x 10^ g 
cm^'^ with a temperature of ~ 3.3 x 10® K. The tremen- 
dous pressure produced by relativistic neutrinos may drive 
materials of such high density to explosion (e.g. Woosley 
& Janka 2005). During the first ~ 10 s of the core col- 
lapse, a power of about ^ 10^^ erg s~^ is released as 
high-energy neutrinos within a radius of ~ 10^ km (e.g. 
Woosley & Janka 2005). The neutrino-electron cross sec- 
tion was estimated to be ~ 10~'*^(E/GeV)cm^ with E 
being the neutrino energy (e.g. Marciano & Parsa 2003). 
During the gravitational core collapse of a SN explosion, 
a typical value of E would be E~20MeV (e.g. Hirata et 
al. 1987; Arcones et al. 2008). Therefore, the neutrino- 
electron cross section is estimated to be ~ 2 x 10"'*'* cm^. 
If we adopt the neutrino luminosity L — 10^^ erg s~*, 
the ratio of neutrino pressure to one electron and the 
iron-core gravity on one silicon nucleus is ^ 10^^. By 
these estimates, neutrino pressure is unable to split the 
silicon layer from the iron core. A pure vacuum void is 
unlikely to appear during the gravitational core collapse 
and the rebound process to initiate a SN. 

During the subsequent dynamic evolution, diffusion 
effects would gradually smooth out any sharp edges and 
the inner rarefied region will be dispersed with diffused 
gaseous materials. Behaviour of this gas may be affected 
by the central neutron star. For example, a relativistic 
pulsar wind is able to power a synchrotron nebula re- 
ferred to as the pulsar wind nebula, which is found within 
shells of supernova remnants (e.g. Gaensler et al. 1999). 
Pulsar winds relate to the magnetic field of pulsars and 
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are not spherically symmetric. However, if the angle be- 
tween magnetic and spin axes of a pulsar is sufficiently 
large and the pulsar is rapidly spinning, the averaged 
pressure caused by a pulsar wind may appear grossly 
spherical. This pulsar wind pressure may also counterbal- 
ance the inward pressure of outer gas and slow down the 
diffusion process. We offer a scenario that a supernova 
remnant makes a transition from a sharp-edged 'void' to 
a quasi-smooth one due to combined effects of diffusion 
and pulsar wind. 

4-3.3 Interstellar Bubbles 

Our ISSV model solutions may also describe large-scale 
nebula evolution around early-type stars or Wolf-Rayet 
stars. Here our overall scenario parallels to the one out- 
lined in Section 4.3.1 for PNe but is different at the am- 
bient medium surrounding the flow systems. For an in- 
terstellar bubble, a central stellar wind collides with the 
ISM (i.e. no longer a dense wind) and gives rise to a 
reverse shock that heats the downstream stellar wind 
zone, and a forward shock that propagates outwards in 
the ISM. Meanwhile, by inferences of radio and opti- 
cal observations for such nebulae (e.g. Carina Nebula 
and Rosette Nebula), the central hot stars arc capable 
of ionizing the entire swept-up shell and thus produce 
huge H II regions surrounding them (e.g. Menon 1962; 
Gouguenheim & BottineUi 1964; Dickcl 1974). The tem- 
perature of a H II region is usually regarded as weakly 
dependent on plasma density, thus an isothermal H II re- 
gion should be a fairly good approximation (e.g. Wilson, 
Rohlfs & Hiittemeister 2008). As the ISM outside the 
forward shock is almost unaffected by the central wind 
zone, we may approximate the static ISM as isothermal 
and the dynamic evolution of gas surrounding interstel- 
lar bubbles can be well characterized by various ISSV 
solutions. Several prior models that describe ISM bub- 
ble shells in adiabatic expansions without gravity, might 
encounter problems. For example, the self-similar solu- 
tion of Weaver et al. (1977) predicts a dense shell with a 
thickness of only ^ 0.14 times the radial distance from 
the central star to the shell boundary, which is indeed a 
very thin shell. However, observations have actually re- 
vealed many ISM bubbles with much thicker shells (e.g. 
Dorland, Montmerle & Doom 1986). The Rosette Neb- 
ula has a H II shell with a thickness of ^ 20 pc; while 
the radius of central void is only 6 pc (e.g. Tsivilev et 
al. 2002). The shell thickness thus accounts for ~ 70% of 
the radius of shell outer boundary, which is much larger 
than the computational result of Weaver et al. (1977). In 
our ISSV solutions, there are more diverse dynamic be- 
haviours of gas shells and outer envelopes. For example in 
Figure 9, four ISSV solutions with static outer envelope 
indicate that the ratio of shell thickness to the radius of 
the forward shock covers a wide range. For these four 
solutions, this ratio is 0.65, 0.98, 0.91 and 0.57 respec- 
tively. 



A type Zj void shock solution with xq = 0.58, ao = 
0.014 and an isothermal shock at Xs — 2.2 may charac- 
terize gross features of Rosette Nebula reasonably well. 
This ISSV solution indicates that when the constant shell 
temperature is 7000 K (somewhat hotter than 6400 K as 
inferred by Tsivilev et al. 2002) and the entire nebula 
system has evolved for ^ 10® years, the central void has 
a radius of ~ 6.0 pc; the forward shock outlines the outer 
shell radius as ~ 22.8 pc; the electron number density in 
the HII shell varies from ^ 8.5 cm^'^ to ^ 12.9 cm^'^; 
the contact discontinuity surface expands at a speed of 
~ 6.1 km s~^ and the forward shock propagates into the 
ISM at a speed of ~ 16.4 km s~^; the surrounding ISM 
remains static at large radii. In the above model calcu- 
lation, an abundance He/H ratio of 0.1/0.9 is adopted. 
Various observations lend supports to our ISSV mod(4 
results. For example, Tsivilev et al. (2002) estimated a 
bit higher shell electron number density as 15.3 cm~^ 
and an average shell expanding velocity of about 8.5 km 
s~^. Dorland et al. (1986) gave an average shell elec- 
tron number density as 11.3 cm~"^. Our calculation also 
gives a shell mass of ~ 1.55 x 10^ Mq, which falls within 
1O,OOOM0 and 16,000Mq as estimated by Menon (1962) 
and Krymkin (1978), respectively. 

To study these inner voids embedded in various gas 
nebulae, diagnostics of X— ray emissions offers a feasible 
means to probe hot winds. The thermal bremsstrahlung 
and line cooling mechanisms can give rise to detectable 
X— ray radiation from optically thin hot gas (e.g. Sarazin 
1986) and the high temperature interaction fronts of stel- 
lar winds with the ISM and inner fast wind with the void 
edge can produce X— ray photons (e.g. Chevalier 1997b). 
We will provide the observational properties of different 
types of voids and present diagnostics to distinguish the 
ISSV types and thus reveal possible mechanisms to gen- 
erate and maintain such voids by observational inferences 
in a companion paper (Lou & Zhai 2009 in preparation). 



5 Conclusions 

We have explored self-similar hydrodynamics of an isother- 
mal self-gravitating gas with spherical symmetry and 
shown various void solutions without or with shocks. 

We first obtain type X ISSV solutions without shocks 
outside central voids in Section 3.3. Based on differ- 
ent behaviours of eigen-derivatives across the SCL, type 
X void solutions are further divided into two subtypes: 
types Xi and Xu ISSV solutions. All type X ISSV so- 
lutions are characterized by central voids surrounded 
by very dense shells. Both types Xi and Xu ISSV so- 
lutions allow envelope outflows but only type Xj ISSV 
solutions can have outer envelopes in contraction or ac- 
cretion flows. 

We then consider self-similar outgoing shocks in gas 
envelopes surrounding central voids in Section 3.4. Type 
Zi void solutions are referred to as the equi-temperature 
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shock solutions with a constant gas temperature across 
a shock front (this is an ideahzation; see Spitzer 1978). 
We also investigate various cases of shocks in type Zu 
void solutions (always a higher downstream temperature 
for the increase of specific entropy). In Section 3.4, we 
have developed the 'phase net' matching procedure to 
search for type Z ISSV solutions with static, expanding, 
contracting and accreting outer envelopes. ISSV solu- 
tions with quasi-smooth edges exist only when gas flows 
outwards outside the shock; all other types of voids are 
surrounded by fairly dense shells or envelopes. 

We have systematically examined voids with sharp 
or quasi-smooth edges for various ISSV solutions. There 
must be some energetic processes such as supernova ex- 
plosions or powerful stellar winds (including magnetized 
relativistic pulsar winds) that account for the appear- 
ance of sharp-edge voids. The denser the nebular shell is, 
the more difficult the void formation is. In other words, 
voids of quasi-smooth edges might be easier to form in 
the sense of a less stringent requirement for the initial 
energy input. In Section 4.1, we point out that the gas 
self-gravity can influence the void evolution significantly, 
especially for the property of regions near void edge as 
shown in Fig. 1121 With the same boundary condition 
of outer medium, ISSV models with and without self- 
gravity can lead to different types of voids, e.g. quasi- 
smooth edge or sharp-edge voids. We suggest that these 
two types of voids may have different mechanisms to gen- 
erate and sustain. Thus the inclusion of gas self-gravity is 
both physically realistic and essential. Besides, we show 
that all voids with quasi-smooth edges are type Z voids, 
that is, with shocks surrounding voids. This indicates 
that shock and expanding outer envelope may well im- 
ply a likely presence of a central void. In fact, observa- 
tions on hot gas flows in clusters of galaxies might also 
be relevant in this regard. For example, McNamara et 
al. (2006) reported giant cavities and shock fronts in a 
distant (z = 0.22) cluster of galaxies caused by an inter- 
action between a synchrotron radio source and the hot 
gas around. Such giant X-ray cavities were reported to 
be left behind large-scale shocks in the galaxy cluster 
MS0735. 6-1-7421. 

Another point to note is that ISSV solutions we have 
constructed are physically plausible with special care 
taken for the expanding void boundary a^o . Void bound- 
ary xo involves density and velocity jumps not in the 
sense of a shock; local diffusion processes should hap- 
pen to smooth out such jumps in a non-self-similar man- 
ner. Nonlinear ODEs ([8]) and ([9]) are valid in intervals 
(0, Xq) and {xq , -l-oo). We have indicated this property 
in Section 3.2 when introducing the concept of ISSV and 
discuss this issue in Section 4.2. There, several plausi- 
ble mechanisms are noted such as powerful stellar winds 
and energetic neutrino driven supernova explosion. More 
specifically, we apply the ISSV solutions to grossly spher- 
ical planetary nebulae, supernova explosions or super- 
nova remnants and interstellar bubbles. 



Our model for planetary nebulae involve three char- 
acteristic interfaces: reverse shock, contact discontinuity 
surface, and forward shock. Steady inner stellar winds of 
different speeds blow on both sides of the inner reverse 
shock and a contact discontinuity surface confines the 
slower downstream wind zone outside the inner reverse 
shock. This reverse shock may be stationary or moving 
(either inwards or outwards) in the laboratory framework 
of reference. The contact discontinuity surface between 
the steady downstream wind zone (on the downstream 
side of the inner reverse shock) and the outer expanding 
gas shell moves outwards at a constant radial speed. Be- 
haviours of outer shocked gas shell outside the contact 
discontinuity are described by type Z ISSV shock so- 
lutions with quasi-smooth edges. Stellar core collapses 
prior to supernova explosions lead to neutrino bursts 
during a short period of time, which might momentarily 
stand against the inward pressure force across the 'void' 
edge and give rise to a sharp-edge 'void' structure. In the 
long-term evolution after the escape of neutrinos, diffu- 
sion effect and outer forward shocks will dominate the 
behaviours of supernova remnants and the sharp edge 
will be smoothed out eventually. In other situations when 
central magnetized relativistic pulsar winds begin to re- 
sist diffusion effects, a quasi-smooth void with shocked 
shell (i.e. type Z ISSV) might also form. 

Similar to PNe, interstellar bubbles may originate 
from strong stellar winds of early-type stars on larger 
scales. We invoke type Z ISSV solutions with quasi- 
smooth void edge to describe the structure and evolution 
of dense shell and outer ISM envelope. In our model, the 
hot shocked stellar wind zone, which is located between 
the reverse shock and the contact discontinuity surface, 
is filled with steady shocked wind plasma. In Weaver et 
al. (1977), the standard Spitzer conduction was included 
to study the shocked stellar wind and shell gas that dif- 
fuse into the shocked stellar wind region. However, the 
stellar magnetic field is predominantly transverse to the 
radial direction at large radii, which will suppress the 
thermal conduction through the hot interstellar bubble 
gas (e.g.. Chevalier & Imamura 1983). As a weak mag- 
netic field can drastically reduce this thermal conduction 
coefficient (e.g. Narayan & Medvedev 2001; Malyshkin 
2001) and a weak magnetic field has little effects on be- 
haviours of the gas shells and nebulae (e.g. Avedisova 
1972; Falle 1975), we do not include thermal conduc- 
tion effect in our model but present a dynamic evolution 
model for interstellar bubbles in terms of a self-similar 
nebular shell sustained by a central steady stellar wind. 
We would note that the existence of a random magnetic 
field may reduce the density 'wall' around a void and 
make the formation of a void easier (Lou & Hu 2009). 
We do not include magnetic field in our model for sim- 
plicity. 

Recent observations of ultraluminous X-ray sources 
(ULXs) show that ULXs may blow very strong winds or 
jets into the surrounding ISM and generate hot bubbles. 
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For instance, ULX Bubble MH9-11 around Holmberg IX 
X-1 has experienced an average inflating wind/jet power 
of ~ 3 X 10'^^ erg s~^ over an age of ~ 10^ years. The 
shock of the bubble travels outwards at ~ 100 km s~^ at 
radius ~ 100 pc. The particle density around the shock is 
~ 0.3 cm-3 (e.g. PakuU & Grise 2008). Approximately, 
this bubble corresponds to a type Z ISSV solution with 
ao ~ 3.8 X 10~^, xq = 0.4 and Xg — 1, and the temper- 
ature of gas shell is ~ 10'' K. Our ISSV model predicts 
a contact discontinuity surface, or interaction surface of 
ULX wind with the ISM, at a radius ^ 40 pc. 

Finally, to diagnose voids observationally, we would 
suggest among others to detect X— ray emissions from 
hot gas. In a companion paper), we shall adapt our ISSV 
solutions for hot optically thin X— ray gas clouds or neb- 
ulae, where we advance a useful concept of projected self- 
similar X— ray brightness. It is possible to detect ISSVs 
and identify diagnostic features of ISSV types which may 
in turn to reveal clues of ISSV generation mechanisms. 
We will also compare our ISSV model with observational 
data on more specific terms. Moreover, projected self- 
similar X— ray brightness is a general concept which can 
be useful when we explore other self-similar hydrody- 
namic or magnetohydrodynamic processes (e.g. Yu et al. 
2006; Wang & Lou 2008). 
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The fcth-order derivative of equation (A2) by the Leibnitz 
rule reads 

i=0 

where CI stands for k\/[i\{k — and ! is the standard fac- 
torial operation. 

Because a{xo) ~ 0, equation (A3) yields a'{xo) ~ 0, 
a"{xo) = and so forth. That is, a{xo) — and equation 
((9]) determine that the arbitrary order derivative of a{x) at 
xo is zero, which is contrary to presumption (Al) based on 
the former assumption that a{x) can transit from zero to 
nonzero across xq of the ZML. Therefore a cannot transit 
from zero to nonzero across the ZML for an isothermal gas. 
For properties of void boundary in a polytropic gas, the in- 
terested reader is referred to Hu & Lou (2008) and Lou & Hu 
(2008). 



B Proof of inequality a"{xo) < 
Equation ((9} can be written in the form of 

[{x - vf - l]a' = a[a - 2(1 - v/x)]{x ~ v) , (40) 
whose first-order derivative d/dx is simply 

2{x - v){l - v')a' + [{x - vf - l]a" 
= {a [a — 2(1 ~ v/x)] -I- a [a' -I- 2v' /x — 2v/x'^'\}{x — v) 

-f a[a - 2(1 - v/x)] (1 - v) . (41) 

We have v{xo) = xo and thus both v'{xo) = and a'(xo) = 
at Xo by equation (|22p . Therefore at Xq, equation (B2) 
becomes 

a"{xo) = -a{xof < (42) 
for a positive density a{xo) > as a physical requirement. 



A Jump of a from zero to 

nonzero value across the ZML 

Let us assume that the reduced mass density a{x) can transit 
from zero to nonzero across the ZML and xo is the transition 
point. For this, we require v{xo) = xo, a{xo) — and when 
x < Xo, a{x) and v{x) vanish. For an arbitrarily small real 
number e > 0, we also require a{xo -I- e) > for a positive 
mass such that there exists a positive integer n for which 
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We now cast equation ((9)1 in the form of 



da 
dx 



2(x — v)/x\{x — v) 
(x — v)^ — 1 



aJ^{x) 



(37) 



(38) 



where ^(x) = [a — 2(x — v)lx\(x — v)l^x — v)'^ — 1]. 

At X = xo, the denominator of f{x) does not vanish 
because v{xo) ~ xo- So J^{x) is a finite continuous analytic 
function near xo- An arbitrary order derivative of J-{x) at xo 
should be finite as well. 
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